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ABSTRACT 

This  document  is  an  introduction  to  the  analysis  of  time  series  using 
autocorrelations.   The  autoregressive,  the  moving  average,  and  the  mixed 
models  for  stationary  and  nonstationary  time  series  —  are  introduced. 
Several  stages  to  analyse  given  time  series  are  discussed.   These  are: 
identification,  preestimation  of  parameters,  the  maximum  likelihood 
estimation  and  diagnostic  checking. 

Numerical  techniques  are  described  in  detail  together  with  some 
examples. 
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Chapter  1   Introduction 

This  document  has  been  prepared  as  an  introduction  to  the  analysis  of 
time  series  -  or  more  precisely,  estimation  of  parameters  in  discrete  time 
series.   A  time  series  is  a  set  of  observations  generated  sequentially  in 
time.   There  are  continuous  time  series  Z(t)  and  discrete  time  series  {Z  }, 
but  we  are  interested  throughout  this  document  only  in  discrete  time  series 
observed  at  equal  interval.   Data  are  measured  at  i,-.,  T0+h,  TQ+2h,  ••••and 
they  are  denoted  by,  for  example,  (Z..  ,  Z  ,  Z„,  ...}  . 

Analysis  of  time  series  is  particularly  important  in  fields  where 
control  and  forecast  are  involved;  observed  data  are  inevitably  conta- 
minated by  random  noise.      Take  the  stock  market,  for  example.   The 
offer  by  sellers  will  be  affected  by  the  previous  price  and  the  most  recent 
change  of  price.   Furthermore  we  must  assume  a  term  which  is  unpredictable. 
Therefore  the  price  X  can  be  described  as 

Xt  =  aSt_1  +  B(Xt.1  -  Xt_2)  +  at 

This  is  a  very  simple  model  and  is  a  linear  combination  of  past  values 
and  random  term.   It  is  called  auto-regressive  process.   After  defining  a 
model,  we  can  estimate  a  and  8  from  a  given  set  of  observations  {X- ,X_.  .  . .  }. 
Then  we  will  be  able  to  investigate  the  effect  of  the  previous  price  or 
change  of  prices  to  the  current  price  of  the  stock  and  also  be  able  to 
predict  the  future  prices.   If  this  model  does  not  give  a  good  enough  result, 
we  may  have  to  extend  the  model  to  include  X  ,  -  X   „  or  a  ,  and  so  forth. 

The  above  example  illustrates  very  descriptively  the  way  a  time  series 
is  analyzed.   Our  main  concern  in  the  following  chapters  is  "how  to  estimate 
good  parameters  from  finite  number  of  data".   "Goodness"  means  better  fitting 
to  data  with  less  number  of  parameters.   {Z  }  is  assumed  to  have  no  periodic 
or  seasonal  characters  in  this  document.   (The  stock  market  usually  shows  a 
seasonal  change,  for  example.) 


From  now  on  a  time  series  is  denoted  by  {Z  }  if  its  average  is  not  zero 
and  {Z  }  if  it  is  adjusted  so  that  its  average  becomes  zero.   The  residual 

random  term  is  denoted  by  (a  }.   It  is  assumed  that  {a  }  is  not  correlated 

2 
and  is  usually  assumed  to  be  Gaussian  E  {a.  }  =  0,  Var  {at}  =  a 

t  t     a 


There  are  several  types  of  models  conceivable  -  autoregressive  models, 
moving  average  models  and  mixed  models,  each  of  which  will  be  introduced  in 
the  following  sections  of  this  chapter,  [1], 


1.1  Autoregressive  Models  (AR  Model) 

In  this  model,  the  time  series  {Z  }  is  described  by  the  following  formula: 

oo 

z  =  y    <j>  z   .  + 


a 
t 


That  is,  Z  is  a  linear  combination  of  previous  values  of  the  process  itself 
and  current  random  pulse.   If  B  is  defined  as  a  backward  shift  operator 
(i.e.  BZ  =  Z  _  ),  (1)  reduces  to 

2     3 
(l-c|>  B-<f>  B  -<j>  B  - )  Zfc  =   at  or  in  short  $(B)Zt  =  at.     (1) 

where  E  {a.  }  =  0.   E  {a  2}  =  a  2. 
£  t      a 

When  we  build  an  AR  model,  $(B)  is  a  polynomial  of  finite  order,  say  p.   In 
particular,  the  AR  process  of  first  order  and  of  second  order 

\  •  h  K-i +  \ 

K   "  *1  K-l  +  *2  Zt-2  +  at 

are  important  in  practice. 


1.2  Moving  Average  Models  (MA  Model) 

This  model  is  described  by  the  following  formula 

Z  _  = 


=  a   +   7   o   a 
t    t    A   j   t-j 
J=l   J 


=  (1+0  B  +02  B2  +  ...)  afc 


=  0(B)  afc  (2) 


-2- 


This  is  called  "moving  average"  because  Zt  is  a  weighted  sum  of  another  pro- 
cess a  .   Hence,  w 
{at}  is  the  input. 


cess  a  .   Hence,  we  can  regard  (Z  }  as  the  output  of  a  linear  filter  of  which 


a 


t 


0(B) 


Z 


t 


Figure  1.1 


0(B)  is  a  polynomial  of  finite  order  q  when  we  build  a  moving  average 
model.   In  particular,  the  process  of  first  order  and  of  second  order 


\=     at  "  +1  Vl 


Zt=      at  "  *1  \-l  "  *2  \-2 


are  important  in  practice, 


1. 3  Mixed  Models 

As  one  can  see,  the  moving  average  models  and  autoregressive  models 
are  in  close  relation  -  one  can  convert  from  one  model  to  another.   In 
other  words  (1)  can  be  written  as 

Z  =  $    (B)  a    and  same  is  true  with  (2). 

An  autoregressive  model  of  the  finite  order,  however,  becomes  of  infinite 
order  in  the  Moving  Average  Model  and  vice  versa.  As  an  extention  of  two 
methods,  it  is  natural  to  consider  models  as  follows: 

$(B)  Zt  =  6(B)  at  (3) 

This  model   is  of   infinite  order   in  both  MA  model   and  AR  model   and   it   is   to 
be  hoped  that   in  some  cases  we  can  describe  a  model  with  less  number  of 
parameters   than  either  MA  or  AR  models.      Mixed  process  with   $(B)    of  order 
p  and  0(B)   of  order  q   is   denoted  by  ARMA    (p,q).      In  practice,   ARMA    (1,1) 


Z     -   <i>,  Z        -,=a      _     <b      a      , 
t        vl      t-1  t  vl     t-1 


is   important. 


-3- 


1.4  Statlonarity  and  Invertibillty 

Let  us  consider  a  first  order  autoregressive  model 

Z  =  a  +  aZ 
t    t     t-1 

If  |  ot  |  >  1,  Z  starting  from  any  value  will  either  grows  bigger  explosively 
or  oscillates  with  increasing  amplitude  as  time  goes  on.   That  is,  a  does 
not  have  enough  power  to  retain  Z  around  its  average  and  the  process 
becomes  almost  deterministic. 

A  model  of  this  kind  is  called  "non-stationary".   This  case  will  not 
be  discussed  in  this  document. 

On  the  other  hand,  if  |a|  <  1  the  effect  of  previous  values  of  Z  die 
out   with  time,  and  a  can  always  play  an  important  role  in  the  process. 

In  general, an  autoregressive  process 

$(B)  Z  =  a 
is  called  stationary  if  the  infinite  series   $   (B)  is  convergent  on  and 
within  the  unit  circle,  or  equivalently,  if  all  the  zeroes  of  $(B)  are  located 
outside  the  unit  circle.   The  explosive  processes  can  be  stated  as  the  cases 
in  which  some  zeroes  of  $(B)  are  located  within  the  unit  circle. 

Sometimes  $(B)  has  several  zeroes  on  the  unit  circle  but  none  are  within 
it.   It  turns   out  that  the  resulting  models  are  of  great  value  in  represent- 
ing some  nonstationary  series  which  may  nevertheless  exhibit  homogeneous 
behaviour.   In  particular,  although  the  general  level  about  which  fluctuations 
are  occurring  may  be  different  at  different  times,  the  broad  behaviour  of 
the  series,  when  differences  in  level  are  allowed  for  may  be  similar.   Many 
series  actually  encountered  in  industry  or  business  exhibit  somewhat  nonsta- 
tionary behaviour  which  can  be  handled  with  this  type  of  model. 

The  polynomial  of  autoregressive  part,  in  such  a  case,  can  be  written  as 

$(B)  =  <KB)  (l"B)d 
where  <J>(B)  is  a  strictly  stationary  operator. 

Thus  a  general  model  is 

KB)  (1-B)d  Zt  =  0(B)  at  (4) 

or   <|>(B)  W  =  Q(B)  a 

where   W  =  (1-B)d  Z  . 
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Homogeneous  nonstationary  behaviour  can  therefore  be  represented  by  a  model 
which  calls  for  the  d   difference  of  the  process.   In  practice,  d  =  0,  1  or  2. 
This  model  is  denoted  by  ARIMA  (p,d,q)  if  cj>(B)  is  of  order  p  and  q(B)  of  order  q, 
ARIMA  is  an  abbreviation  of  Auto  regressive  integrated  moving  average  model. 
The  process  is 


t    1   t-1  p   t-p    t   1   t-1 


0   a 
Uq   t-q 


with  W  =  V  z  . 


The  general  autoregressive-integrated-moving  average  process  may  be  generated 
from  white  noise  a  by  means  of  three  filtering  operations  as  shown  below, 


<j)_1(B) 

,_d,  -1 

Moving 

average 

filter 

Stationary 
auto regressive 
filter 

KV    ) 

white 

N 

Nons  tat  ionary 

Summation 

filter 

>■ 

noise 

e 

t 

w 
t 

Z 

at 

0(B) 


Figure  1.2 


There  is  a  corresponding  property  to  stationarity  in  the  moving  average  method, 
That  is  called  " invert ibility".   Let  us  consider  a  first  order  moving  average 
model  to  illustrate  the  basic  idea  of  invertibility . 


Zt  =    (1-0B)    afc 


(5) 


Expressing    a  's   in  terms  Z    's   yields 


a     =  (1+0B  +  0  B      +  .    . 

fc  2 

or     Z     =  a     -  0Z  -  0     Z     „ 

t  t  t_i  t-2 


)   z 


(6) 


The  process  (5)  itself  is  definitely  stationary  for  any  0.   However,  if 
|0|  > 1,  the  weights  in  the  corresponding  autoregressive  process  (6)  grow 
bigger  as  j  in  Z  _ .  increases.  We  can  avoid  this  ridiculous  situation  by 

putting  |0|  <  1.   In  general,  a  MA  process  Z  =  0(B)  is  called  "invertible" 

-1  ii 

if  0   (B)  is  convergent  for  |B|  <  1,  or  if  all  the  zeroes  of  0(B)  are  outside 

of  the  unit  circle. 

As  we  shall  discuss  later,  a  convergent  expansion  for  a  is  possible  when 
the  above  conditions  are  not  satisfied,  but  only  in  terms  of  present  and  future 
values  of  the  process.   The  invertibility  is  needed  if  we  are  interested  in 


-5- 


associating  present  events  with  past  happenings. 

In  the  following  discussion,  both  conditions  for  stationarity  and  in- 
vertibility  are  to  be  imposed  on  parameters  of  models. 


-6- 


Chapter  2   Linear  Stationary  Models 

The  first  section  is  devoted  to  the  autocorrelation  function  p,  , 

which  is  a  very  powerful  tool  in  analyzing  the  time  series.   Linear 
stationary  and  non-stationary  models  which  were  introduced  in  Chapter  1 
are  discussed  in  more  detail  in  the  following  sections.   The  formulae 
which  are  derived  in  this  chapter  are  to  be  used  in  the  subsequent 
chapters  in  estimating  the  parameters. 


2.1   Autocorrelation  Function 

When  a  time  series  is  stationary,  that  is,  in  a  state  of  statistical 
equilibrium,  we  can  define  the  mean  value,  variance  and  autocorrelation. 
Since  the  process  is  stationary,  the  distribution  p(Z  )  is  the  same  for  all 
times  t  and  can  be  written  p(Z).   Hence  the  process  has  a  constant  mean 


00 

y   =  E  [Zt]  =^/%(Z)dZ 


and  a  constant  variance 

O    2    =    E  [(Z  -  y)2]  =^(Z-y)2  p(Z)  dZ 
z         t        -°° 

If  we  have  N  observations  {Zn ,  ....Z  }  we  can  estimate  y  by  the  mean 

In 

of  the  time  series 

1   N 

2 
and  the  variance  o"   by  the  simple  variance  of  the  time  series 

N 

%2=Ui(zt-Z)2-  <2> 

The  condition  of  stationarity  also  guarantees  that  the  joint  distribution 

p  (Z  ,  Z   )  is  the  same  for  all  t..  and  t„  if  they  are  constant  apart,  i.e. 
12 

i   i  2 
if  a  time  series  satisfies  the  condition  E  {|Z  |  }  <  °°  for  all  t  and 

E  {Z    Z  }  =  R(t)  does  not  depend  on  s,  it  is  called  stationary  in  wide  sense, 


-7- 


Then,  we  can  define  the  autocovariance  y     at  lag  k  by, 
Y=  cov  [Zt,  Zt+k]  =  E  [ (Zt-y) (Zt+k  -y  )]. 


If  Z   is  an  ARIMA  process,  the  autocovariance  can  be  derived  from  the 
generating  function.   This  is  to  be  discussed  in  the  next  section. 

The  sample  autocovariance  C,   is  the  estimate  to  Yi,: 


Ck  =  S  \   »!-*>  (ZI+k  "S>- 


(3) 


The  autocorrelation  p  at  lag  k  is 


Yk     E  [(Zt-y)  (Zt+k-y)] 


Y0 


(4) 


where 


p0 


=  1. 


Th 


e  corresponding  sample  autocorrelation  is  Yi<-=  C,  /C  , 


Covariance  matrix  Tn  associated  with  a  stationary  process  for  observa- 
tions made  at  n  successive  times  is 


Tn  = 


Yo      Yl      Y2 Yn_! 

2 

-1  P1 

P2.- 

pn-l" 

Yl    Yo                Y    9 

z 

pl  x 

v                     n-2 

?2 

In 

N           YX 

's 

t 

X 

\ 

• 

Yn-lYn-2            Yl   Y0 

_pn-l 

N 

The  matrix   is  non-negative  definite,   but   takes   0   for  trivial   cases. 
(Example)      Calculate  Z,    C   ,    y  ,    of   the  data  given  below. 
Z  =    {47,    64,    23,    71,    38,    64,    55,    41,    59,    48} 

Z  =  Yq      (47  +  64  + )   =   51 

Z  =    (-4,    13,    -28,    20,    -13,    13,    4,    -10,    8,    -3). 

co  =  To    ^2  +  12>2  +  •*•  +  a2)  =  189*6 


(5) 


C,    = 


1  =   10      ((-*)    x  13  +  13x   (-28)   +   +  8  x    (-3))    =  -1497 


=  -A     =   -0.79 
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Estimated  values  of  the  autocorrelations  are  also  stochastic  variables 
and  their  standard  errors  are  very  important  when  we  want  to  know  whether  or 
not  the  theoretical  autocorrelations  are  truncated  at  some  lag  q  and 
larger.   Naturally  the  error  becomes  larger  if  size  of  the  data  is  not 
sufficiently  large.   Bartlett  [  2]  gave  a  formula  for  the  variance  of 
estimated  autocorrelation  when  it  is  known  that  only  p. , p„ . . . . p  are  non 
zero.   That  is 

1       q    2 
Var  [yk]  -  -  {1  +  2^   p  u}    if  k  >  q  (6) 

u=l 

Therefore,  we  can  substitute  y      ....  y     for  pn  ....  p  ,  respectively  to 

1       q       1       q 

see  whether  the  assumption  that  p,  =0  (k  >  q)  is  valid  or  not.   This  method 
will  appear  in  the  identification  procedure  of  parameter  estimation. 


2.2  Autocovariance  Generating  Functions  of  ARIMA  Processes 

There  is  a  very  convenient  formula  to  give  any  order  of  autocovariance 
function.   Suppose  an  ARIMA  process  is  given  as 

$(B)  Zfc  =  0(B)  at. 

Operation   $      (B)    from  the   left  yields 

Zt    =    *_1(B)0(B)at    =    *  (B)    at 


=      I  *.    Bj    a     =      I        *, 


a 


j-0        J  £      j=0      J     t-j 


00  00 


^-EtZtW  ■  J0    ^'i'j  E[Vi  \W 

a2  6 
a        i,    j-k 

00  »        «J 

=   "a"      Z   n   Tj    fJ+k-  <7> 

J=o 

00 

2  2  7  2 

Y     =   a       *  a        =     L        y.      .  (8) 

o  z  a  .  j 

j=0 


Especially 
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The    autocovariance  generating   function  is  defined  as 

00  00  00 

Y(B)    =   X      Y     Bk  =  o'2y        I        V.    ¥...     Bk 

{  kLcoYk  a     £.*  3      3+k 


T 
=0    for-    all     v    <    n 


Since  ¥1  =0   for  all  k  <  0, 


00  00 


Y(B)   -  a  2     I       I  Vh  Bh"j  <h  =  J+k> 

a     j=o  h=0  J 

00  00 

j=0     J  h=0  h 


=  a   2   $_1(F)    0(F)    •    S>   -1    (B)0    (B)  (9) 

a 


(Example)      MA(1)    Process 

Zt  =    at  -    9at_1  =    (1-0B)    afc 

Y(B)=    (1-0B_1) (1-0B) 

=   1+  02    -   0B"1   -   0B 

,o=l+02 

•Vk=   0  (k  >   2). 

2.3     More  About  Autoregressive  Models    (Stationary) 

An  autoregressive  model  of  order  p,    AR   (p) ,    is  denoted  by 

K  '   *1   \-l+  *2   Zt-2  + +  *p  Vp  +    \  (10) 


or   in  short      $(B)    Z     = 


\ 


2 
where     $(B)   =   1   -   $     B  -    <j)     B     - 
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2.3.1  Autocorrelations  of  Autoregressive  Models  (Yule-Walker  equation) 

Multiplying  (10)  by  Zt_k  and  taking  the  expectation  at  each  term  yields 


Vt,  =  *i  Yv_,  +  <f>9  Y,  9  +...  +  *  Y,  _ 


1  'k-1   Y2  'k-2 


P  'k-P 


k  >  0 


(11) 


because  E  [Z  k  at]  =  0  for  k  >  0.   Note  that  the  expectation  vanishes  when  k  <  0, 
Hence 

=  A    ^       +  A    ^       -4-  An 

(12) 


pk  =  *1  pk-l  +  *2  pk-2  + 


P  k"P 


As  can  be  seen  from  (12)  the  autocorrelations  for  AR  processes  is  infinite  in 
extent.   If  we  take  the  first  p  equations  in  (12)  and  write  them  in  matrix 
notation,  we  get: 

1   P- 


'1   p2 

pl   X   P2 
1 


p-1 


p-1 

p-2 


Ca     "1 

+1 

*? 

*p 

(13) 


(13)  is  usually  called  Yule-Walker  equation  and  is  very  useful  in  the  preliminary 

estimation  of  <J>'s  from  sample  autocorrelations. 

2 
The  variance  a   is  obtained  from  (10)  by  taking  its  expectation  after  multiplying 
z  2 

it  by  Z  .   Since  Z   is  related  to  a  .  E  {Z  a  }  =  a 
t  t  t      t   t      a 


Hence 


a 
a 


1-  *x  Pr  *2  P2- 


.-<j>  P 
P  P 


2.3.2  Examples 

(1)   Zt  =  izt_±  +   at 


(AR(1)  process) 


This  process  is  called  Markov  process. 

P,_  =  + 

P,    =  *k        (k  >   1) 

R  2  2 

2_  a a 

°z    "  !"    pl*       =      1  -*2 

|<j)J  <     1    (stationarity   condition) 


(14) 
(15) 

(16) 
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(2)      Zfc  =    <|>1  Zt_1  +   <f>2   Zfc_2  +   at  (AR   (2)   Process) 


The  stationarity  condition  yields 


♦x  +  *2  <  1 

4>2  -  *]_  <  1 
-1<  (j>2  <  1 

Pk  "  *1  Pk-1  +  *2  Pk-2 


Po 

Pi 

= 

1 

l-*2 

p2 

= 

*2 

<J>     and   <J>?   are  estimated   from  Yule-Walker  equations: 


P2  -    Pl2 

*2   =      1  -    pj 


2.4  The  Partial  Autocorrelation  Function 

There  is  another  powerful  function  which  helps  us  check  whether  a  process 
is  autoregressive  or  not.   Suppose  a  time  series  is  described  by  AR  (p)  model. 
Then  the  autocorrelation  function  p.,  although  infinitely  extended,  is  a 
linear  combination  of  p  non  zero  functions  of  autocorrelations.   For  example, 
in  AR  (1)  process  Z  =  <|>Z.  •,  +  a,.   (See  2.3.2)   p  =  <j>  is  the  only  independent 
function.   So  let  us  investigate  how  we  can  exploit  this  fact  to  identify 
the  autoregressive  models. 

Let  us  assume  that  the  process  is  of  the  order  £  so  that 

\"  ♦uit«  +**2V2  + +  t»Vi+at  (17) 

From    (12)   we  obtain 

pi    =      *Upj-l+  hi  Pj-2+ +  *U  PJ-»  (18) 

j   =  1,2...* 
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where  we  can  solve  <f>„  ,'s  in  terms  p.  's  by  the  Yule-Walker  equations. 

We  solve  these  equations  for  £  =  1,2,3,....  successively  to  get  a  sequence 
(  <f>   ,   <J>99,   <J)„„,....).   Since  <j>   is  the  coefficient  of  the  last  term 
in  (17),  it  has  a  cut  off  after  lag  p  if  the  original  process  is  of  degree 
p.   Quenouille  showed  [10]  that  the  variance  of  the  estimated  <j>,  ,  from  n 
data  if  k  >  p  is  given  by 

var  [<f>kk]~  ~    for  k  >  p  +  1. 


Hence  we  can  know  the  degree  of  autoregressive  models  by  plotting 


Jkk 


against   k   and  by   checking    when  it    starts   to   be   submerged  within   the   region 
defined  by    Jvax    [<k,]    =  n  ~2*      An   example   is   given  in  Fig.    2.1. 


A 


-1 
n   2 

0 

-1 

n    "2" 


Fkk 


m 


a 


Figure  2.1 

An  example  of  typical  behavior 

of  the  partial  autocorrelations 


There  is  a  very  efficient  algorithm  to  compute  the  partial  autocorrelatio 
function  without  solving  Yule  Walker  equations  at  each  step.   This  is  a 
recursive  method  in  which  the  estimate  for  AR(p)  can  be  computed  from 
those  of  AR  (p-1)  and  autocorrelations.   The  general  recursive  formulae  are 


Vi, 

j 

*Pj  "  Vi 

,  p+1 

*P: 

p-j+1 

p+1 

p 

I 

i=l 

Vi 

rp+i-.i 

-i-1 

*    • 
PJ 

r  . 
3 

j  =  1,2,3, .. .p 


(18) 


(19) 


These  formulae  are  actually  used  in  the  program  USID. 


-13- 


(20) 


2.5  More  About  Moving  Average  Models 

A  moving  average  model  MA(q)  is  denoted  by 

Z  =  a  -  0  a    -  0       0  a 

t    t    1   t-1    2   t-2  q  t-q 

=  (1-  0,  B   -  0  Bq)  a 

1  q      t 

-   0(B)  at 

2.5.1  Autocorrelation  function  of  Moving  Average  Models 
The  autocovariance  function  of  a  Mflfc  (q)  process  is 

Tk  -  I  [  (»t  -  \   at_1  0q  St_q)  (at_k-Sl  at.k_x . . .  -Sqat_k_q) ,   (n) 

Since        E[at_katHl]    -o^   ^    ,    (21)   yields 

Hence       Yo  -  (l+e  2  +  e2  + +e  2)    a  2 

i  z  q 

and  y,    =        (6.    1^    9...    +9.   9,  ,„  + +  0  0  .    -   2..  O0. 

k  k        1     k+1  2     k+2  q-k     q)    u      (k=l,2 q)    (22) 


Y.    -  0  (k  >   q.) 


We  see  that   the  autocorrelation  function  of  a  MA(q)    process  has   a   cut   off 
at    lag   q.      The  difficulty  with  moving  average  models   is   that   the  equations 
(22)    are  not  linear  if  we  are  to  solve  them  with  regard  to   0's.      Therefore, 
they  must  be  solved  by   some   iterative  methods.      The  algorithm  will  be 
found   in  3.3  and  the  actual   program  in  Appendix. 

2.5.2     Examples 

(1)      Z     =  a      -   0     a.„   .  (MA   (1)    process) 

t  t  i     ^-1 

invert ibility   condition     -1    <   0  <  1 

2  2 

Y     =    (1+  0..    )    a 
o  la 


k  >    2 


pl 

= 

-01 

1+012 

pk 

= 

0 

e1  =      (  -i  +  'i-4 Plz  )/2P;L 


*      =    -  0k  {      1  -  el2 } 

9kk  1  -,         ©  2(k-l) 

1 


(2)   Z  =  a  -  0  &  -  e  a  (MA  (2)  process) 

t     t     1   t— ±     Z   t— z 


Invert ibility 


02  +  0X  <  1 
02  -  0X  <  1 


-l  <  e2  <  i 


2     2     2 
Y  =  (1+0,  +  e9  )  o 
o       12a 


-0!     (1-    0?) 

Ml " 

2            2 
1+  01     +  02 

"    09 

p2  - 

2           2 
1  +   0X   +  02 

pk  = 

0                      k  >    . 

2.6  Duality 

One  will  notice  the  close  relation  between  AR  models  and  MA  models. 
As  has  been  discussed  in  Chapter  1,  an  AR  model  of  finite  order  becomes  a 
MA  model  of  infinite  order  and  vice  versa.   Therefore,  the  autocorrelations 
of  a  MA  (12)  process  have  a  cut  off  at  lag  k,  whereas  its  partial  auto- 
correlations are  extended  to  infinity  because  the  process  is  equivalent  to 
an  AR  process  of  infinite  order.   Note  also  that  the  dual  of  the  stationary 
condition  in  AR  models  is  the  invertibility  condition  in  MA  models.   The 
duality  will  become  more  obvious  in  Mixed  Models. 

2.7  More  About  Mixed  Models 

We  have  noted  in  Chapter  1  that  in  some  cases  mixed  models  require 
less  number  of  parameters  than  MA  models  or  AR  models  because  the  mixed 
models  have  the  characteristics  of  both. 

We  employ  the  mixed  models  as  follows: 

Z  =  0n  Z  .  + +  0  Z    +  a-0n  a    -  0   a     .0  a        (23) 

t     1   t-1  p   t-p     t   1   t-1     2   t-2     q   t-q      v   ' 


or  in  short 


*(B)  Zt  =   0(B)  afc  (24) 


This  process  is  referred  to  as  ARMA  (p,q) 
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(24)  defines  a  stationary  process,  provided  that  all  zeroes  of  $(B)  are 
outside  the  unit  circle,  and  defines  an  invertible  process  if  all  zeroes 
of  0(B)  are  outside  the  unit  circle.  We  assume  that  these  conditions  are 
satisfied.   Remember  that  $(B)  has  some  zeroes  on  the  unit  circle  when  we 
refer  to  ARIMA  processes. 

2.7.1  Autocorrelation  function  of  mixed  models 


On  multiplying  (2.7.1)  by  Z    and  taking  the  expectation  yields 

Yk  =  *1  Vl  +  ••••  +  *p  Vp+  Yza  (k)"ei  YZa(k-l)-.--0q  Yza(k-^ 
where  y      (k)  =  E  [Z  ,  a  ]  is  the  cross  covariance  function  between  z  and  a. 
Obviously  y      (k)  =  0  if  k  >  0. 

Z  Si 

Therefore 

Y,  ■  6,  y,  .  +  e.  Yi  o+ +  ©  Yi        Ck  >  q+l) 

1  k    1  'k-1    2   k-2        p  'k-p  n 

Hence 

pk=  *1  pk-l  +  *2   pk-2+--"+  ^p  Pk-p 


Or 


♦(B)  P,  =  0  (k  >  q+l) 


Thus,  autoregressive  parameters  are  separated  from  moving  average  parameters 
if  we  take  lag  k  >  q+l.   We  can  derive  a  modified  MA  (q)  process 

~Zt=\-hK-l-    *2^t-2------Vt-q 

if  we  know  all  the  autoregressive  parameters.   This  fact  is  used  in  the 
preliminary  estimation  of  parameters.  (Chapter  3). 

2.7.2   (Example) 

ARMA   (1,1)   process 

Practically,    this   process   is   very   important. 

Z     -   cL    Z     ..    =    a    -    0.      a- 
t  1      t-1  t  1        t-1. 


or 


(1-   ^B)    Zfc      =      (1-G.jB)    at    . 
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Autocorrelations  are 

Yo  -  ^  Yjl  +  Oa2   -   01  yza  (-1) 

Yl  =  h   Yo  "  91  aa2 


Yk  =  *!  Yk_!  k  >_2. 


It  is  easy  to  find  v   ("!)• 

za 

Y   (-1)  =  (^  -  0-)  a„2 
za  1     1    a 


Therefore  „ 


Yo    = 

(1-(|)101)(<J>1-   91) 

2 

aa 

2 

'1  - 

-*x2 

aa 

(l-cj)101)((j)1  -  01) 
Pl  =  2 

i+e1    -  2^  ex 

P2  =   Vl 


2 
from  which  <))..  ,  <j>9,  and  a   can  be  calculated, 


Finally,  it  is  worthwhile  to  summarize  the  typical  behavior  of  autocorrela- 
tions and  partial  autocorrelations  of  three  models.      If  the  process  is 
AR  then  a.c.  has  exponential  decay  or  damped  oscillations  whereas  p.a.c.  has 
a  sharp  cut  off.   If  the  process  is  MA,  then  the  situation  is  reversed.   In 
MIX  processes  both  a.c.  and  p.a.c.  decay  exponentially  or  have  damped  oscillation 
(there  are  four  combinations  depending  on  the  value  of  parameters).   (See 
Section  3.1  for  a  summary.) 
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Appendix  Al.   Backward  Representations  of  Time  Series 

Suppose  that  W,  is  represented  by  the  linear  model 

(KB)  Wfc   =  0(B)  at  , 

where  the  zeroes  of  <J>(B)  and  0(B)  lie  outside  the  unit  circle.   It  is  very 

interesting  to  note  that  the  autocovariance  generating  function  y(B)  (§2.2), 

given  by, 

y(B)  =  $(F)    0(F)   $   (B)  0(B) 

is  invariant  under  the  exchanges  B  -*■  F  and  F  ->  B.   Thus  the  stochastic  process 

<KF)  Wt  =  0(F)et  (1) 

(even  though  its  meaning  is  not  clear  as  yet)  has  the  identical  covariance 
structure  as  the  original  stochastic  process 

KB)  Wt  =  0(B)  afc.  (2) 

e   in  (1)  is  a  sequence  of  independently  distributed  random  variables  having 

2     2 
mean  zeros  and  variances  o       =    a 

e     a 


The  stochastic  process   (1),  means  that  W  can  be  expressed  entirely  in  terms 

of  future  W's  and  e's  and  is  a  stationary,  invert ible  representation.   We  refer 

to  it  as  the  backward  form  of  the  process.  We  will  use  it  in  the  estimation  of 

the  parameters. 
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.Appendix  A2 .   Spectral  Properties 

The  autocorrelation  function  is  not  the  only  tool  to  analyze  a  time 
series  with;  sometimes  the  power  spectrum  is  very  useful.   Actually  in  the 
fields  such  as  oceanography  or  communication  engineering,  the  power  spectrum 
is  used  rather  than  the  autocorrelation  function. 

In  this  document,  however,  the  power  spectrum  appears  only  in  Chapter  6 
to  check  the  periodicity  of  a  time  series,  therefore  a  very  brief  review  is 
presented  here.   The  reader  should  refer  to  [3]  for  more  details  on  this  subject 

1)   Periodgram 

In  the  analysis  of  the  spectral  properties  we  regard  a  given  time 
series  as  composed  of  sine  and  cosine  waves  of  different  frequencies.   First 
we  define  periodogram  of  a  time  series.   That  is,  if  the  number  of  the 
observation  is  N  =  2q+l.  (odd) 

Z  =  a  +  )    (a.  cos  2ir  f.t  +  b.  sin  2tt  f.t)  +  e 
t    q   .  ,    i-        i     i        it 
i=l 

where  f .  =  —  is  the  ith  harmonic  of  the  fundamental  frequency  1/N.   Fourier 

coefficients  a  ,  a.,  b.,  are  estimated  as 
o    1    x 


O      L 


2   N  N 

ai  =  N  I      Zt  cos  2tt  f.t,  b.  =-  I      Z   sin  f.t     t-1,2,3, . . .N. 

t=l  t=l 

N-l 
The  periodogram  then  consists  of  q  =  — ^—  values 

N    2     2 

i±  =  2  <ai   +  bi  >         ±  =  i.2.»-.q 

When  N  is  even,  we  set  N  =  2q  and  obtain 
a  =  Z 

0  2  ? 

a.  =  —  )   Z  cos  2tt  f.t       i=l,2,...  q-1 

1  N  f\   t         l  »  »     -l 

N 
1 
a 


q  ■  s  lx  ^  \ 

2   N 
bi  =  N  I        Zt  S±n  2U 


t=l 


i 


b   =  0 

q 
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N        2 
I.  =  f  <ai  +bi  >      i=l,2,...q-l 

1=1  (.5)  =  Na  2 

q         q 

Note  that  the  highest  frequency  is  .5  per  time  interval. 

I.  indicates  the  relative  intensity  of  the  ith  frequency  component.   If 

Z  contains  a.cos2fr  f.t  or  b.  sin  2tt  f  .t  or  both,  the  periodogram  shows  a 
t  ill         i  r      ° 

peak  at  i.   On  the  other  hand,  if  Z  =  a  +  e  ,  the  periodogram  is  constant. 
This  is  why  we  can  discuss  some  periodic  features  of  a  time  series  by  using 
its  spectral  properties. 

2)   Sample  spectrum  and  power  spectrum 

I.  assumes  a  discrete  argument.   We  define  1(f)  to  have  a  continuous 

1  i 

argument  f  and  to  take  values  I.  at  f .  =  —  .   This  new  function  is  called 

the  sample  spectrum.   There  is  a  very  important  relation  between  1(f)  and 

the  autocovariance  function  C.  of  Z  . 

k  t 

Theorem.      The  sample  spectrum  is   the  Fourier  cosine  transform  of  the  auto- 
covariance  function: 


N-l 

1(f)   =   2    {Cn  +  2      I     C     cos   2tt  f,  }  0   <   f   <  T 

K=l      k  k  Z 


Proof: 


1(f)   =  |  (af2+  bf2)   =  |  (af+ibf)    (af-ibf) 

N      * 
=  2  df  df 

where  N 

d-  =    a-ib_  =  -     I      (Z      cos   2rr   ft   -  iZ     sin  2irft) 
r  r        r        JM      -..        t  t 

N 
=  |  ^    (Zt-i)   e"2*  flt 

•••     1(f)   -  |      I  I      Vt-V    »t-  -  «    e  -*»«<t-0 

t  =  l  ^ 

1      X  t'=l 

0      N-l      N-k  _       ... 

=  1      I        I       (z-  -z")    (Z-+w  "  Z>  e 
NK=(n-l)j=i  J  J+k 

N-l  „    ... 

-   2  I  Ck  e   ~2-fk 

k=-(N-l)       K 

N-l 
=   2      {C  +2        J        C,     cos    2tt    fk}    . 


k-1 
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I[f]  is  an  estimate  which  was  obtained  from  N  samples  {Z. ,Z_,....Z  }. 

LI  n 

It  is  a  stochastic  variable.   We  take  the  expectations  value  of  1(f), 

N-l 
E  {Iff]}  =  2  {E[CJ  +2    T   E  [C.  ]  cos  2Trfk} 
U       k=l      R 

As  N  goes  to  infinity,  E  [C .  ]  tends  to  the  theoretical  autocovariance  function 
y..   Therefore,  we  can  define1 a  power  spectrum  p(f)  as  a  limit  of  E  {l[f]}  : 

p(f)  =  lim    E  {1(f)}  =  2{y  +2Z  y     cos  2n   fk}  . 

O        K. 

n  ->  °° 
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Appendix  A3.   Conditional  Expectation 

The  conditional  expectation  of  Z  at  time  k,E  [Z  ]  is  the  expectation 
given  complete  historical  knowledge  up  to,  but  not  beyond  time  k.   For  the 
simplicity  in  notation,  let  square  brackets  imply  that  the  conditional 
expectation  at  time  t  is  to  be  taken.   For  example, 

Ut+J,l  =  E  [a£+4] 

Some  properties  are: 

(1)  [Z  _.]  =  E  [Z   ]  =  Z  j  =  0,1,2,.... 

U   J        *.  J  J 

That  is  to  say,  standing  at  time  t  the  expected  values  of  Z's  that  have 
happened  already  are  the  values  they  have  actually  realized. 

(2)  [Z  ,.]  =  prediction  of  Z,..  at  time  t. 

(4)   [at+j]  =  E  [at+j]  =  0. 

i.e., at  time  t,  the  expected  values  of  the  a's  that  have  yet  to 
happen  are  zero. 
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Chapter  3   Identification  Of  Models 

Now  we  are  prepared  to  tackle  practical  problems;  given  a  time  series, 
we  can  estimate  the  parameters.   There  are  three  stages  for  analyzing  the  time 
series.   Firstly  we  must  identify  the  type  and  the  order  of  the  process, 
we  can  obtain  crude  approximation  to  the  parameters  at  this  stage.   Secondly 
we  must  improve  the  estimate  of  the  parameters.   And  finally  we  must  check 
whether  or  not  the  estimated  model  really  fits  the  given  time  series. 
This  chapter  is  concerned  with  the  first  step,  i.e.  identification  of  the 
model.    Our  principal  tools  are  the  autocorrelations  and  the  partial  auto- 
correlations. 

It  is  inevitable  that  identification  and  estimation  overlap.   To  identify 
a  model,  we  must  make  some  inference  on  parameters.   But  in  this  stage  the 
estimate  is  inexact.   Our  aim  here  is  to  identify  the  model  to  a  given  time 
series  in  the  following  way: 

1.  to  difference  Z  as  many  times  as  is  needed  to  produce  station- 

arity,  reducing  the  process  under  study  to  the  ARMA  process 

<|>(B)  W„  =  0  +  ©(B)  a 
to         t 

where  1J  =  (1-B)d  Z„ . 
t  t 

2.  to  identify  the  resulting  ARMA  process. 

Programs  are  also  given  in  appendices.   They  are  explained  in  the  text. 

3.1  Use  of  Autocorrelations  and  Partial  Autocorrelations  in  Identification 
Procedure 

3.1.1.  Identifying  the  degree  of  differencing. 

We  have  seen  in  Chapter  2  that  the  autocorrelations  die  out 
quickly  if  none  of  roots  <f>(B)  =  0  are  close  to  the  unit  circle.   Therefore, 
if  the  plotted  autocorrelations  do  not  show  quick  decay,  it  indicates  that 
further  differencing  is  necessary.   In  practice  the  order  of  differencing  is 
0,1  or  2  and  it  is  usually  sufficient  to  inspect  the  first  20  or  so  estimated 
autocorrelations. 

3.1.2.  Identification  of  resultant  stationary  ARMA  process 

Having  tentatively  decided  the  order  of  differencing,  we  next 
study  the  general  appearance  of  the  estimated  autocorrelations  and  partial 
autocorrelations  of  the  appropriately  differenced  series.   The  characteristic 
behavior  of  these  functions  summarized  in  Chapter  2  is  very  helpful  to 
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identify  the  type  of  the  model.   Both  autocorrelations  and  partial  autocorre- 
lations should  be  plotted  for  first  20  lags  or  so.   This  can  be  nicely  done 
by  line  printers  or  CRT  displays.   (See  the  appendix). 

Some  important  characteristics  of  the  models  are  summarized  in 
Table  3.1. 

TABLE  3.1 


order 


(1,  d,  0) 


(o,  d,  1) 


\k 


estimates 


decays  exponentially 
only  (()-.  1  non-zero 


*1  =  P1 


only  p.,  non-zero 

exponential  dominates 
decay 

-Ql 

pi  "  i+e1  2 


admissible  region 
(stationarity, 
invert ibility) 


-1  <  ^  <  1 


-1  <  0  <  1 


order 


(2,  d,  0) 


(0,  d,  2) 


mixture  of  exponential 
or  damped  sine  wave 


only  p   and  P2 


non-zero 


'kk 


estimates 


only  <}>   and  4>__  non-zero 


♦x- 


pi  (1-pi' 
i-P,2 


dominated  by  mixture 
of  exponential  or 
damped  sine  wave 

-e1  d-e2) 

1  =     2   2~ 
1+0,  +02 


* 


P2  -Pl 


2    .     2 
1-  PX 


-02 


2   2 
1+01  +02 


admissible  region 


•1   <   <|>2   <  1 


-1   <  ©2  <  1 


«j>1  +<J>2  <  1 


02  +  01  <  1 


<j>2  "   +!   <  1 


e2  -  ex  <  i 


-24- 


rkk 


estimates 


admissabl«: 
region 


(1,  d,  1) 
decays  exponentially  from  first  lag 
dominated  by  exponential  decay  from  first  lag 


Pi"" 


(i-o1<l>1)(*1  -e1) 


1+0, 


P2  -  Px  *x 


-1  <  ^  <  1, 


2  4>1  e1 


-1  <  0  <  1 


3.1.3.   Standard  error  of  the  autocorrelations  and  the  partial  autocorre- 
lations 

Usually  the  estimated  autocorrelations  and  partial  autocorrelations 
fairly  fluctuate  around  the  axis  and  take  moderately  large  values  after  they 
once  decay.   Therefore  we  must  have  a  broad  view  in  identifying  models.   Some- 
times we  may  have  multiple  models  until  estimation  stage  or  even  until  diag- 
nosis stage.   For  large  lags  there  are  some  formulae  which  give  clues  to  the 
standard  errors  of  those  functions.   For  example,  if  the  theoretical  auto- 
correlations are  0  after  lag  q,  the  standard  error  of  the  estimated  autocorre- 
lations is  given  as 


ff[rk]  -  -  (1+2  (r±2   +  r/  + 


1/2 


+  rq» 


k  >  q   (see  (2.1)). 


As  for  the  partial  autocorrelation  function,  the  estimated  function  has  the 
standard  error 


0  [  V^     k  >  p 

if  the  theoretical  function  is  0  for  k  >  p. 


(2.4) 


Therefore  we  can  draw  la  or  2a  lines  in  the  graph  of  those  functions  to 
indicate  whether  they  are  negligible  or  not. 
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3.2        Description  of  program  USID 

This  program  is  for  the  identification  of  the  models. 

1)  Input  parameters 

Z:   time  series  to  be  analyzed  (saved) 

N:  number  of  data 

CD:  maximum  lag  of  autocorrelation 

L:  maximum  lag  of  partial  autocorrelations 

D:   order  of  higher  difference 

2)  Output 

MEAN:   mean  value  of  Z  (Z) 

VAR:  variance  of  Z  (C  ) 

ATCR:   Autocorrelations 

o[y   ]:      Standard  error  of  autocorrelations 

cr[cf> ,,  ]:  Standard  error  of  partial  autocorrelations 

3)  Subroutine 

CORRECOPRINT :   print  out  the  function 

4)  Function 
Read  as  data. 

Take  the  higher  difference  as  specified  (D) 
Calculate  the  mean  value,  variance,  auto  and  partial  auto- 
correlations. 

Calculate  the  standard  error  of  them. 
Plot  the  result,  via  a  line  printer 

3. 3  Preliminary  Estimates  of  Parameters 

Now  that  we  have  an  appropriately  differenced  time  series  W  =  V  Z  , 
we  can  make  initial  estimates  of  parameters.   Formulae  for  lower  models  are 
given  in  Table  3.1. 

There  are  general  formulae  which  work  for  any  order  and  which  are  very 
suitable  for  computer  programming.   They  will  be  shown  in  the  subsequent 
sections. 

3.3.1  AR  Parameters 

Autoregressive  parameters  are  obtained  from  Yule-Walker  equations: 

<f>  =  R  ^r   .  (2.13) 

P   P 
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This  is  a  system  of  linear  equations  and  is  solved  by  any  standard 
method.   (It  is  recommended  to  have  a  subprogram  to  solve  a  system  of  linear 
equations  in  one's  programs.)   It  can  be  proved  that  the  autoregressive 
parameters  obtained  from  Yule-Walker  equations  for  the  estimated  autocorre- 
lations approximate  the  fully  efficient  maximum  likelihood  estimates.  (Proof  in  [1].) 


The  initial  estimate  of  the  residual  variance  is 
2 


ca  =  Co  (1-  ^   r±   -  ^   r2 


P   P 


3. 3. 2  MA  Parameters 

In  general  the  equations  for  moving  average  parameters  are  non-linear 
and  much  harder  to  solve  than  AR  parameters.   More  explicitly,  they  are 


-G,  +6.  9,  ,.  + 0  .  0 

r       k   1   k+1 q-k  q 

k=       ,  .  :  2  .  :  2  .  .  : 2 


k  =  1,2, 


1  +  9   +  02  +  . 


+  6 


CQ    =  aa      (1+  Q±   +....+  0q  ) 

This  system  of  nonlinear  equations  is  solved  iteratively.   The  method 

discussed  here  is  a  Newton-Raphson's  algorithm  which  converges  quadrat ically, 

Suppose  we  have  the  estimated  autoc ©variances ,  C~,  C.  ,  C-....C  . 

0   1   2     q 

Procedures : 


1'      T0  =/  C0' 


2.   Solve 


T-.  =  ....=  T   =0 

i        q 


i+l 


i+l 


x  i+l 

q 


t  i 

q 


t  _1  f1 
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where 


T  = 


V  Tr 


q-l 


T0  Tl 


Tl   V 


f  = 


q-j 


•  -    I      t. 


j-o 


T.,  .   "  C. 

1+J    J 


for  x    iteratively  until  | f , |  's  become  sufficiently  small, 


3.   Finally 


"Tl/T0 
"T2/T0 


_"VToj 

'  2     2 
aa  =  TQ 

This  algorithm  has  been  given  by  Wilson,  [4]. 

3. 3. 3  Mixed  Process 

It  is  often  found  that  W  =V  Z   is  most  economically  represented  by  a 
mixed  ARMA  process 

<KB)  Wt  =  9(B)  At. 

The  procedure  to  estimate  parameters  of  mixed  ARMA  process  is  a  mixture  of  the 
two  described  in  3.3.1  and  3.3.2. 
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Suppose  we  are  given  the  first  p+q+1  autocovariances  C.  ,C   . . . . ,  C    , 


Then  we  have  three  stages 

:egr 
We  discussed  in  2.5.1  that  if  the  lag  is  greater  than 


(1)   The  autoregressive  parameters  <J>-  <f>  are  estimated  from 


q-p+1'  q+p 

q  then  only  autoregressive  parameters  are  related  to  autocovariance  functions. 

Therefore,  we  can  solve  Yule-Walker  equation  for  <f>.'s. 

(2)  Using  the  estimates  <J>.'s  obtained  in  (1),  the  first  q+1  auto- 
covariances C/  (i  =  0,1.... q)  of  the  derived  series  W  "  =  W  -  A,  W  ,  ,..-d>   W 

1  t     t    1  t-1     p   t-t 


We  do  not  have  to  recalculate  W  's  to  obtain  C's.   The 


are  computed. 

following  formula  yields  directly  the  derived  autocovariances. 


Cj  "  I     J     *i*k  C|j+i-k|  •    ♦o'"1'  (J  "<>.•.«)■ 

J    i-0  k=0         IJ     ' 

(3)   Finally  the  derived  autocovariances  C^  C'        are  used  in  an 

o       q 

iterative  calculation  to  compute  initial  estimates  of  the  moving  average 

"  2 
parameters  0, ,0„,..,0  and  the  residual  variance  a 
L      2  a  a 


3.4  Description  of  Program  USPE 

This  program  is  for  the  initial  estimates  of  parameters 
1)   Input  parameters 


W 
N 
P 

Q 

CD 

MEAN 
VAR 
ATCR 
EPSIL0N 

2)   Output 

PHI 

THETA 

SIGMA 


time  series 

number  of  data 

order  of  AR  process 

order  of  MA  process 

maximum  lag  for  autocorrelation  function 

mean  value  of  w  (W) 

variance  of  W  (G_) 

autocorrelation 

criterion  for  the  convergence  of  Newton-Raphson' s 
method 

AR  parameters 
MA  parameters 
Residual  Variance 


-29- 


3)  Subroutine 

GAUSSLU     Gaussian  Elimination 

4)  Function 

Compute  AR  parameters  by  Gauss  elimination. 

Compute  derived  autocovariances. 

Compute  T  matrix  and  f . 

Newton-Raphson's  algorithm. 

Compute  MA  parameters  and  residual  variance, 

Print  out. 
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Chapter  4   Maximum  Likelihood  Estimates 

The  identification  procedure  has  led  to  a  tentative  formulation  for 
the  model.  We  now  need  to  obtain  efficient  estimates  of  the  parameters. 
The  method  of  the  maximum  likelihood  function  (essentially  the  method  of 
least  squares)  is  very  suitable  for  this  purpose. 

The  first  section  is  an  introduction  to  the  maximum  likelihood  func- 
tions.  Readers  familiar  with  the  subject  can  skip  this  section.  Other  sections 
of  this  chapter  are  devoted  to  the  usage  of  the  maximum  likelihood  functions 
to  estimate  parameters. 

4 . 1  The  likelihood  functions 

The  method  of  maximum  likelihood  is  often  a  very  suitable  means  of 
obtaining  the  required  estimates.   In  many  cases  the  errors  of  the  esti- 
mates can  be  smaller  than  those  which  would  arise  from  any  other  method 
of  estimation.   The  concept  of  maximum  likelihood  estimators  is  introduced 
by  some  examples. 

(Example  1  -  discrete)   Two  urns  are  available,  each  containing  black  and  white 
balls.  Urn  A  has  3  black  balls  for  every  white  ball,  and  urn  B  has  3  white 
balls  for  every  black  ball.   One  of  the  urns  is  chosen  at  random,  and  then 
three  balls  are  picked  up  with  replacement  from  the  urn.   The  probability 
of  getting  r  black  balls  is 

P3(r)  =  (3r)    pr  (l-p)3"r  (1) 

where  p  =  3/4  for  urn  A 
1/4  for  urn  B. 

So  far  the  problem  is  that  of  probability.   Now,  we  want  to  estimate  p  on 
the  basis  of  the  observed  example.   (Note  that  the  problem  goes  backward!) 
We  can  evaluate  (1)  for  all  possible  values  of  r  and  for  the  two  values  of 
p.   That  is, 


r 

0 

1 

2 

3 

P(rl3/4) 
P(r|l/4) 

1764 
27/64 

9/64 
27/64 

27/64 
9/64 

27/64 
I/ 64 

(*)  This  section  is  based  on  "Lectures  on  Elementary  Statistics  and 
Probability"  by  D.  J.  Hudson,  CERN«  Report  #63-29. 
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Suppose  we  in  fact  observe  no  black  balls,  i.e.  r=0.  We  then  argue  that 
the  observed  data  are  27  times  more  likely  to  have  occurred  if  p  =  1/4 
than  if  p  =  3/4.  We  may  conclude  that  we  must  have  chosen  the  urn  B. 
We  would  say  that  the  most  likely  estimate  of  p  =  1/4  if  r=0  or  1,  and 
3/4  if  r=  2  or  3. 


We  can  extend  the  problem  to  the  case  where  there  are  9  urns  each 
of  which  have  the  required  mixture  of  black  and  white  balls.   (Such  as 

Pl  = 


,  f 2  -  u,^. 

*q 

-  \j.  j ) . 

iae  Lduit 

:  uecuuu 

r 

0 

1 

2 

3 

P(r  |1/10) 

.729 

.243 

.027 

.001 

P(r|2/10) 

.512 

.384 

.096 

.008 

P(r  |3/10) 

.343 

.441 

.189 

.027 

P(r  |4/10) 

.216 

.432 

.288 

.064 

P(r  |5/10) 

.125 

.375 

.375 

.125 

P(r  |6/10) 

.064 

.288 

.432 

.216 

P(r  |7/10) 

.027 

.189 

.441 

.343 

P(r  |8/10) 

.008 

.096 

.384 

.512 

P(r  [9/10) 

.001 

.027 

.243 

.729 

For  a  given  possible  value  of  r,  we  choose  p  which  is  most  likely  to  have 
given  rise  to  data.   (Underlined  ones). 

(Example  2  Continuous) 

Suppose  that  the  random  variable  Z  has  a  continuous  distribution  with 
density  function, 
f(z|  0) 
where  0  is  a  single  parameter  which  we  want  to  estimate.   The  distribution 
of  a  sample  of  size  n  is 


p  (z. ,  z_ z  )  ndz.  -  n   f  (z.  le)  az.. 

12      n      i      ,      i  '     1 

Once  the  Z  's  have  been  observed,  the  likelihood  of  0  is  defined  by 
in 


(2) 


£(  o|z)  -  n  f  (z  |e  ) 
i=i    1 


(3) 


which   is  of  the  same  form  as  P  but  in  which  Z's  are  fixed  and  0  becomes 
variable. 
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The  log  likelihood  of  6  is  defined  by 

L  (0|Z)  =  Zlog  f  (ZJ0).  (4) 

The  estimated  0  is  the  solution  of  j^     =  0  or  —  =  0.   The  reason  why  L 

is  used  is  that  the  exponential  form  appears  very  often  in  statistics. 

If  the  normal  distribution  with  unknown  variance  and  mean  is  assumed,  i.e., 

-(Z-y)2  /    2 


f(Z|v,a  )  =  y^-2  e  2o  (5) 


2 


A(y 


,a2|z)  =  (a2)-?exp-  I    (Vy)  /2a2  (6) 


1=1         2 
n   (Z.-  y)Z 


L(y,c/|Z)  --f  log  aZ  -  £  — i_~  (7) 

1  i=l  2a  z 


|t-0    +  ^  I(Z.-y)  =  0   .'.  y=Z 


(8) 


2 

3L    n       n    X  -i.  1        \fv        ^2    n    •    ~2    ^  "  ZJ) 
-r-2=  0   ->  -  -r-   -r  + — r      )(z.  y)   =0   .  .   a   =  (9) 

8  or  2    a2  ^q  1  n 

In  many  examples,  for  moderate  and  large  samples  the  log-likelihood 
function  will  be  unimodal  and  can  be  adequately  approximated  over  a  suffi- 
ciently extensive  region  near  the  maximum  by  a  quadratic  function.   The 
second  derivatives  in  such  cases  provide  measures  of  "spread"  at  its 
maximum  of  the  likelihood  function  and  can  be  used  to  calculate  approximate 
standard  deviations  for  the  estimates. 

4. 2  Use  of  the  Maximum  Likelihood  Functions  to  Estimate  Parameters  of  Time 
Series 

4.2.1 

Now  let  us  discuss  how  to  use  the  maximum  likelihood  functions  to 
estimate  the  parameters  in  ARIMA  models.   The  idea  is  to  calculate  the 
residual  terms  by  using  the  pre-estimated  parameters  and  to  use  the  maximum 
likelihood  function  for  Gaussian  distribution  with  pre-estimated  variance* 
Suppose  we  have  n  +  d  observations.   The  ARIMA  model  (p,d,q)  can  be  written 
as 

t    t    1   t-1      p   t-p    1   t-1    2  T:-2      q   t-q 


where    W  =  VdZ   w  =  W.  -  E[W  J 


(d  >  0  -*  E[W  ]  =  0). 


(10) 
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It  can  be  shown  (proof  omitted,  See  [  l])  that  the  log-likelihood 
function  is  given  by 

L  (4,6,0  )  -  f  (*,6)  -  nln  a   -   S(*te?  UU 


2a2 
a 


2 


S(4>,0)  =  I        [a J  0,  <)>,W]   is  called  the  unconditional  sum  of  squares 

t=-oo 

function.   [a  |<}>,0,W]  denotes  the  expectation  of  a  conditional  on  <f>,0  and  W. 
(See  Appendix  3  of  Chapter  2.) 

f(<}>,0)  is  important  only  for  small  n  and  L  is  dominated  by  the  last 
term. 

Therefore,  the  parameter  estimates  obtained  by  minimizing  the  sum  of 
squares,  least  square  estimates  usually  provide  very  good  approximation  to 
the  maximum  likelihood  estimates.   The  unconditional  sum  of  squares  are 
calculated  by  computing  [a  ]'s  recursively. 

4.2.2   General  procedure  for  the  unconditional  sum  of  squares 

Suppose  that  the  W.'s  are  generated  by  the  stationary  forward  model 
cj)(B)  W  =  0(B)  a.  (12) 

whereV  Z  =  W  and  W  =  W  -u  (u=0) .   The  corresponding  backward  formula  is 

<KF)  Wt  =  0(F)  afc.  (13) 

1)     From   (12)   and  (13)  we  obtain  the  dual  set  of  equations 

♦(B)  [Wt]  =  0(B)  [at]  (12-) 

*(F)  [Wt]  =  0(F)  [et]  (13-) 


We 


first  compute  [e  ]'s  by  (13').   Several  e  's  beyond  e   are  assumed  to 
L  t  n 


be  zero  to  start  off  the  sequence.   [e  ],  [e.  ,] [e.  ]  are  thus  obtained. 

n     n-1       1 

[e  ],  [e   ]....  are  zero  because  they  are  distributed  independent  of  W. 

2)   We  then  calculate  [W_,]  by  using  [W,]'s  and  [e  ]'s.   Since  there  are 
autoregressive  operators  (Jk's  [W_.]  are  extended  to  infinity  although  [e_,]'s 
vanish.   But  the  condition  of  stationarity  ensures  that[W  ,]'s  decrease  and 
becomes  essentially  zero  beyond  some  point. 

.  3)  After  [W  . ]'s  vanish,  we  can  compute  the  forward  sequence  [a, ]'s  by 
(12").  [a_.]'s  =  0  beyond  the  point  where  [W_,]'s  -  0.   Thus  [a,]'s  are 
obtained  beyond  n. 
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v    2 
4)   Compute  S  (<J>,0)  =  Zfat]  • 

An  illustrative  example  is  given  below  to  understand  the  general 

procedure.  In  table  (4.1) [W  ]  is  the  given  data.   An  assumed  model  is 

(1-0. 3B)  W  =  (1-0-7)  a 
Then,  (12^ )and(13^ )  can  be  written  as 

[at]  =  [Wfc]  -  0.3  [Wt_1]  +  0.7  [a.+1],  and 

[e  ]  =  [W  ]  -  0.3  [W  x]  +0.7  [et+1]  respectively. 


where  [W  ]  =  W  , 


(t  =  1,2, 


n) 


TABLE 

4.1 

t         [ 

*ti 

0.7    [at_ 

.1]-0.3[Wt_1] 

[wt] 

-0.3[Wt+1] 

0.7fet+1] 

[et] 

-4 

A 

■  0  ! 

!o 

-3 

-0.01 

0.00 

0.00 

-0.01 

(2)      0.01 

!  ° ' 

,' o 

-2 

-0.04 

-0.01 

0.00 

-0.03 

0.03 

!  ° ' 

1  0 

-1 
0 

-0.11 
-0.36 
-1.20 

-0.03 
-0.08 
-0.25 

0.01 
0.03 
0.09 

-0.09 
-0.31 
-1.04 

0.09 
0.31 
0.60 

'0J 

1.64 

1  0 
'0 

1      .    ■ 

1 

1.47 

-0.84 

0.31 

"2.0 

-0.24 

0.58       ' 

N  2.34 

2 

1.23 

1.03 

-0.60 

i    0-8 

0.09 

-0.06 

0.83 

3       (3) 

0.32 

0.86 

-0.24 

1-0.3 

0.09 

0.13 

-0.08 

4 

0.02 

0.23 

0.09 

-0.3 

,             0.57 

-0.09 

0.18 

5 

-1.80 

0.01 

0.09 

,-1.9 

-0.09 

1.86 

-0.13 

6 

-0.39 

-1.26 

0.57 

'   0.3 

-0.96 

3.32 

2.66 

7 

2.84 

-0.27 

-0.09 

1    3.2 

i           -0.48 

2.02 

4.74 

8 

2.63 

1.99 

-0.96 

1.6 

0.21 

1.08 

2.89 

9 

0.66 

1.84 

-0.48 

;-0.7 

-0.90 

3.14 

1.54 

10 

3.67 

0.46 

0.21 

1    3.0 

-1.29 

2.78 

4.49 

11 

5.97 

2.57 

-0.90 

,'4.3 

-0.33 

;  o ; 

3.37 

12 

N 

3.99 

4.18 

-1.29 

1    1.1 

i I 

(1) 

; 

j   init: 

.ally   given 

1  0 

or  as 

sumed 

S(0.3,  0.7)  =  I  [a  J   =  89.2 

t— 4 

We  can  perform  the  second  iterative  cycle,  but  in  most  of  the  cases  it  is  not 
necessary;  the  convergence  is  very  rapid.   In  this  example,  the  second 

r     2 

iteration  yields   S  =•     2,[et]      =   89-3. 


-35- 


4. 3  Graphical  Study  of  the  sum  of  squares  functions 

It  is  very  helpful  to  use  some  graphical  displays  to  determine 
the  region  of  the  global  minimum  of  the  sum  of  squares  functions. 


S(0)< 


Smin 


1  -  parameter  case 


2  -  parameter  case 


If  S  is  a  function  of  two  parameters,  a  contour  map  or  a  "cross  section  graph" 
are  very  useful  to  understand  the  characteristics  of  S. 

Some  Remarks 

(1)  We  must  examine  the  total  characteristics  of  the  maximum  likelihood 
functions.   We  should  note  that  there  are  more  than  one  minimum  in  some  cases. 

(2)  It  is  possible  to  detect  overfitted  case  by  checking  whether  the 
minimum  is  located  on  the  boundary  of  the  parameter  space. 

(3)  Sometimes  there  are  constraints  on  the  domain  of  parameter  space. 
In  such  cases,  the  maximum  likelihood  function  must  be  minimized  within  the 
limited  domain. 


4.4  Variances,  Covariances  and  Confidence  Region 

Near  the  minimum  point,  the  maximum  likelihood  function  can  be  expanded 

by  Taylor's  theorem: 

k        k 
L(?)   =  KS.o2)-   L    (a.o2)   +  y     I       I     £      (3.-3.) (3.-3.)   + (14) 

3.  3.  "i=l      "f  =1  J         J 


where   I. .   =    , 
ij        93 


32L 


i   J 


,    k  =  p+q    (the  number  of  parameters) 
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The  variance  matrix  is   then  defined  as: 

V(3)    =      {'E[£..]}      _1   =    2   o    2    {S      }   _1  (15) 

ij  a  ij 

where  2    0,QiTT>. 

S  •  •    "     I    a      II     <k  x  k  ^trix) 

i j  3   3.93. 


(It   is   assumed  that   the  size  of  samples  are  moderate  and   large  so  that   S   is 
the   dominant  part  of  L   in    (14) .      It   can  be  shown  that 

;   2   a   SJC61  (  f  omitted,    see    [1]).  (16) 

an 

Furthermore,    the   confidence  region  can  be  determined   for  a  given   confidence 
level   e   in  the   following  way: 

-   I  E[Z..]    (3.-  3.)    (3.    -   3.)<      Xc  2    (k)  (17) 

2        2 
where  x   (k)  is  X  -  distribution  function  with  degree  of  freedom  k.   Using 

the  above  relation,  the  confidence  region  for  a  given  confidence  level  e  is 

x/(k) 
S(3)  <  S(3)    (1  +  — )  •  (18) 

4.5  On  Parameter  Redundancy 

Theoretically,  the  model 
<j>(B)  Wfc  =  0(B)  at 
is  identical  with  the  model 

(1  -  aB)  <f>(B)  W  =  (1-aB)  0(B)  a  . 

But  in  the  estimation  procedure  this  redundancy  causes  serious  difficulties. 
Therefore  we  should  avoid  the  situation  where  redundancy  or  near-redundancy 
occurs.   Let  us  take  an  example  in  ARMA  (2,1)  model: 

(1-1.3B+0.4B2)  W  =  (1-.5B)  a  . 

or  (1-.5B)(1-.8B)  W  =  (1-.5B)  a  . 

which  is  identical  with 
(1-.8B)  W  =  a  . 

A  near  redundant  case  can  be,  for  example 

(1-.4B)  (1-.8B)  W„  =  (1-.5B)  a  . 
t  t 
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There  are  more  than  one  combination  of  parameters  which  yield  similar 
[a  ]'s  and  consequently  similar  likelihoods  can  be  found;  a  change  of  para- 
meter value  on  the  left  hand  side  can  be  nearly  compensated  by  a  suitable 
change  on  the  right  hand  side  to  give  the  same  likelihoods. 

The  sum  of  squares  contour  surfaces  in  the  3-dimensional  parameter 
space  will  be  cylinders  rather  than  ellipsoids  and  a  line  of  near  least 
squares  solutions  will  be  obtained  rather  than  a  clearly  defined  point  minimum. 

The  identification  procedure  in  Chapter  3  is  very  important  in  the  sense 
that  it  enables  us  to  avoid  such  disastrous  situations.  Here  is  an  example  in 
which  the  direct  cancellation  of  factors  occur  in  AKMA  (1,1)  process: 

(1-4>B)  Wt  =  (1-0B)  at. 

If  <}>=  0  then  W  =  a  , 

that  is,  W  becomes  a  white  noise  process  itself.   The  sum  of  squares 
function  is  constant  on  0=  <J>.   In  practice  the  identification  technique 
will  reveal  any  such  situation  easily.   Note  the  following  things: 

(1)  We  should  avoid  mixed  processes  containing  near  common  factors  and 
we  should  be  alert  to  the  difficulties  that  can  result. 

(2)  We  will  automatically  avoid  such  processes  if  we  use  identification 
and  estimation  procedures  intelligently. 

4.6  Description  of  Program  MLES 

Input  parameters   W  ,  (time  series) 


Procedure 
Output 


2 
°a   ,  0,  <l>  (pre-estimated  parameters) 

RESIDUAL  (residual  terms) 
S  (least  square  function) 


This  program  computes  the  sum  of  the  squares  function  from  given  data  using 

the  pre-estimated  parameters.   First  the  residual  terms  a  are  reproduced 

2 
and  then  the  least  squares  function  is  computed  by  summing  up  a   's. 
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Chapter  5 Nonlinear  Estimation 

In  the  method  described  in  Chapter  4,  we  must  compute  the  sum  of  the 
squares  function  S(3)  many  times  in  the  parameter  space  around  the  pre- 
estimated  parameters  3nJ  to  locate  the  point  which  gives  the  minimum  of  S(3),  i.e. 
S(3).   This  is  very  time  consuming  when  there  are  many  parameters  to  be 
estimated,  although  it  gives  global  behavior  of  the  function  S(3)  and  is  less 
dangerous.   We  can  improve  this  by  using  some  nonlinear  methods.   The 
classical  Gauss-Newton's  method  is  introduced  in  Section  1  and  its  modifi- 
cation by  Marquardt  is  discussed  in  Section  2.  Computation  of  the  derivatives 
of  S(3)  is  discussed  in  Section  3.   Comparison  of  two  methods  is  also  given. 

5.1   Linearization  of  the  models  (Gauss-Newton's  method) 


Let  us  take  an  ARMA  (p,q)  model 

(J>(B)  w  =  0(B)  a  t  =  0,1,2,. ,.,n. 

Let  3  denote  9's  and  <J)'s.   It  is  a  (p  +  q)  dimensional  vector.   (Let  us  put 

k   =  p  +  q.)   We  can  regard  [a  ]'s  as  continuous  functions  of  3  and  expand 

them  in  Taylor's  series  up  to  the  first  order  around  a  given  set  of  initial 
values  3n» 

[at]  -  [a  ]  +  J     (B.  -  S.0)  Jt  +  0  <(B.  -  6,/)    (1) 
U     1=1  l 


where 


[a  ]  =  [a J  w,  e  ]. 
c0 


In  matrix  notation,  (1)  is  written  as 
[t]    =  -  X  (3  -  3Q)  +  [aQ] 

where 

+  i 

[a]    and    [a„]    are  N   dimensional  vectors    and   x..    =   -  — — 

1J  J 

Then  S(3)    is    computed   as    follows: 


.->- 


S(3) 

t=l 


■     1     at    (3)2   -    [tf    .    [a] 


=    [lQ]t.     [aQ]    -   2    [t0f   X   (3   -   3Q)   +   (3   -   3Q)    XC   X(3   -   3Q) 
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go 

By  putting  -tt-  =  0,  we  obtain  a  system  of  equations 

Xfc  X(3   -   3Q)    =  X   [aQ]  (2) 

or 

N  k  N 

t=0     j=l  c'x     C»J        J  Ju  t=0        ,J 

N  N 


By 


putting  P..   -      I     XX  and  g     -     J     X         a         we   can  write   (2)    as 

XJ         t=0  t=0 


P    (63)   =  g  (3) 

where  63  =  3  -  3n  is  the  correction  to  the  parameters.   Therefore  we  can 
obtain  the  improved  values  for  3  by  solving  (3)  with  regard  to  63. 

(3)  is  usually  solved  after  scaling,  [5].  That  is, 


P*  =  (P..)  = 


JO 


/p  .  .   /p.  . 


*        §-• 


g*  =  (g,)  = 


Then  a  new  system  of  equations 

P   63  =  g  .  (3') 

is  solved  exactly  (by  Gauss  elimination).     63   in  (3)  is  obtained  from  63 

-1 

\      ->  * 

63   =/  \     63 

^22 


* 


or  63,  =  63,  / 


i     "i  '^7 
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Since  the  Taylor's  expansion  (1)  is  correct  only  near  3~  (it  gives  only  the 
tangential  plane  of  S(3)  at  3=  3n)  ,  a  single  adjustment  will  not  immediately 
produce  the  least  square  values;  several  iterations  are  necessary.   Naturally 
convergence  is  faster  if  reasonably  good  initial  guesses  are  used.   This  is 
why  we  have  been  bothered  so  much  in  making  good  preliminary  estimation  of 
parameters  in  3.3. 

5. 2  A  Modification  by  Marquardt 

Sometimes  Gauss-Newton's  method  turns  out  to  be  very  unstable.   Mar- 
quardt improved  its  convergence  property  by  increasing  all  the  eigenvalues 
of  P*  in  (3')  by  a  fixed  positive  amount  [6].   The  theory  of  his  method  is 
not  shown  here,  but  the  algorithm  is  briefly  discussed. 

After  obtaining  the  scaled  system  (3*) ,  the  matrix  P*  is  changed  so 
that 

P„  =  P*  +  XI 

M 

where  X  is  a  positive  number.   Then 

■>         -*  (4) 

PM  63*  -  g*  K   J 

M 

is  solved.   It  can  be  proved  that  the  new  parameters  3  =  3n  +  63  will  lead 
to  a  new  sum  of  squares  S(3),  which  is  always  smaller  than  S(3)  for  a 
sufficiently  large  X,  unless  3^is  already  at  a  minimum  of  S(3).   Our 
strategy  to  find  X  is  as  follows: 


(1)  If  S(3")  <  S(3),  the  parameter  correction  63   is  tested. 
If  max  1 63.  I  <  e   where  e  is  some  prescribed  small  number  such  as  10   , 
convergence  is  assumed.   Otherwise,  3   is  modified  and  X  is  multiplied  by 
t(t  <  1)  and  computation  continues  with  new  parameters.   (i.e.  new  P  and  g 
are  calculated). 

(2)  If  S(S')  >  S(3),  X  is  divided  by  x  (note  that  X/x  >  X  ) 

and  computation  is  continued  to  solve  the  system  of  linear  equation  (4)  with 
new  X  until  reduced  sum  of  squares  S(3)  is  found  or  X  becomes  unreasonably 
large.   It  is  suggested  by  Marquardt  [6]  and  Bard  [7]  to  put  X  =  . 01  initially. 
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5. 3  Numerical  estimates  of  the  derivatives 

It  is  necessary  to  calculate  the  derivatives  X   .  to  construct  the 

i»  J 
system  of  equations  (2).   There  are  two  methods;  the  direct  method  and  the 

numerical  method.   The  former  is  to  make  use  of  some  recursive  formulae  and 

will  not  be  discussed  in  this  document.   The  latter  is  essentially  the  same 

as  deriving  residual  terms  in  4.2.2.   We  calculate  both 

[a  |  ex,  62,  ...  efc]  and  [a  |   3-^  $2,...B  +60  ,  Bfc] 

and  then  obtain  X^  ,  by 
t,j 

[aj  3X,  e2, B±  +  66± efc]  -  [at  |ei,  62 Bfc] 


t,j  60± 

This  method  has  the  following  advantages: 

(1)  applicability   is  universal. 

(2)  the  same  subprogram   (RESIDUAL)    as    in  4.2.2   can  be  used. 


5.4     Comparison  of  the  two  methods 

Both  methods   have  been  implemented   and  tested.      With   the  test   example 
given  in  Appendix    (197  data)    Gauss-Newton's  method  turned   out   to  be 
unstable;    oscillation  was   observed.      Marquardt's  method,    on  the  other 
hand,   was   stable   and   converged  very   quickly.       (Table  5-1). 
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TABLE  5-1 


197  Data  (See  the  appendix) 


Values  of  S(3) 

"^^^rnodel 
Iteration 

(1,  0,  1) 

(0,  1, 

1) 

N  -  G 

Marq. 

N  -  G 

Marq. 

1 

19.264 

19.260 

322.07 

320.76 

2 
3 

19.255 
19.254 

19.254 

318.35 
318.36 

318.35 
318.32 

conv. 

4 
5 

19.254 
19.254 

318.36 
318.37 

318.31 

conv. 

6 
7 

19.254 

318.37 
318.37 

conv. 

8 

318.38 

9 

318.38 

10 

318.39 

11 

318.38 

12 

318.40 

13 

318.39 

14 

318.40 

15 

318.39 

16 

318.41 

17 

318.39 

18 

318.41 

19 

318.39 

20 

318.41 

Note:   Convergence  is  checked  by  max  |<5$.|,  not  by  6S(3). 
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5.5  Description  of  program  MLQT 


Input 


W 

N 

e 
p 
Q 

mean 

£ 

6 

I  max 

T 

X 

H   max 


time  series 
size  of  w 


preestimated  order  and  parameters 


0) 

criterion  for  convergence  of 
used  to  calculate  the  derivatives 
maximum  number  of  iterations 


) 


explained  in  5.2 
a  bound  for  X 


Output 


S 

e 
x 


sum  of  the  squares 
parameters  estimated 

variance  covariance  matrix  (to  be  added) 


Subprograms 


GAUSSLU 
RESIDUAL 


Functions 


1. 
2. 

3. 


Computes  derivatives  using  RESIDUAL. 

Computes  improved  values  of  parameters  using  GAUSSLU. 
(Marquardt's  algorithm.) 

2 
Computes  S  =  Ea   and  X  (to  be  added) 
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Chapter  6   Model  Diagnosis 

Having  estimated  the  parameters,  we  must  examine  whether  or  not  the  model 
which  we  have  built  really  fits  the  given  data.   One  way  of  checking  the 
models  is  called  a  method  of  overfitting,  that  is,  to  estimate  the  parameters 
in  a  model  somewhat  more  general  than  that  which  we  believe  to  be  true. 
This  method  assumes  our  knowledge  on  the  direction  in  which  the  model  is 
likely  to  be  inadequate.   Therefore  it  is  not  general  and  will  not  be  dis- 
cussed here.    The  second  way  of  checking  the  models  is  to  make  use  of  the 
white  noise  term    which  we  have  computed  to  obtain  the  sum  of  the  squares 
function  S;  to  check  whether  a  s  are  sampled  from  Gaussian  population. 
There  are  two  methods;  that  is,  the  one  which  employes  (1)  the  autocorrela- 
tion functions  and  (2)  the  cumulative  periodogram. 

6.1   Autocorrelation  Check 

Suppose  a  model 

d)(B)  W  =  0(B)  a   with  W  =V  Z  has  been  identified  and  the  maximum 
T     t         t        t     t 

likelihood  estimates  (  <j>,  0)  obtained  for  the  parameters. 

As  in  Chapter  5,  the  quantities   a  =  0   (B)   <)>(B)  W  .    t=l,2,...,N 

2 
are  supposed  to  have  Gaussian  distribution  with  mean  0  and  the  variance  a 

i  a 

Therefore  it  is  to  be  expected  that  study  of  the  a  s  indicates  the  adequacy 

of  the  model. 

2 
It  is  known  in  X  -  theory  that  the  quantity 

l(\   ~    at)2/  a*  (1) 

A      1     A  9  9 

where  a  =  —  I   a  is  distributed  as  x  (n-1),  if  a  's  are  sampled  from  N(0,a   ). 

Therefore,  if  n  is  sufficiently  large,  then  the  quantity  (1)  should  be 

2 
approximately  N-1  and  the  ^   testing  with  some  confidence  level  is  available. 

In  practice,  however,  this  criterion  is  too  crude. 

Hence,  more  sophisticated  diagnoses  are  necessary.   How  about  using  the 
autocorrelations  r(k)  obtained  from  a  's?   If  ar's   are  white  noise,  the  auto- 
correlations r(k)  k  >  1  will  be  normally  distributed  around  0  with  standard 
deviation  ~—  .    Unfortunately,  this  diagnosis  is  not  accurate  enough  in 
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practice,  either;  unless  the  lag  k  is  sufficiently  large,  this  criterion  is 
underestimated  and  dangerous  to  use. 

The  best  diagnosis  to  exploit  the  autocorrelation  functions  is  called 
"portmanteau"  test  in  which  r,  (a)'s  are  taken  as  a  whole  to  test  inadequacy 
of  the  model  rather  than  individually.   Suppose  that  we  have  the  first  K 
(say,  20)  autocorrelations  r,  (a)  from  an  ARIMA  (p,d,q),  then  it  is  possible 
to  show  [8]  that,  if  a  's  are  sampled  from  Gaussian  population  the  quantity 

Q  =  n  I       xl    (a) 

k-1   k 

is  approximately  distributed  as  X  (k-p-q),  where  n=  N-d  is  the  number  of 

2 
W's  used  to  fit  the  model.   Hence,  we  can  examine  x  -testing  with  regard  to 

Q. 

6. 2  Cumulative  Periodogram  Check 

This  is  the  only  section  in  this  document  where  the  power  spectrum 
(periodogram)  is  used  in  practice.   We  have  assumed  in  the  beginning  that 
there  is  no  seasonal  or  periodic  nature  in  the  data.   But  in  actual  cases 
there  will  be  many  data  which  involve  some  periodic  characteristics.   There- 
fore we  discuss  in  this  section  a  way  of  detecting  periodicity  in  a's.   The 
autocorrelation  function  does  not  indicate  sensitively  such  periodic  nature. 
The  power  spectrum,  on  the  other  hand , is  very  suitable  to  detect  periodic 
patterns  buried  in  white  noise. 

The  periodogram  is  defined  in  A2  of  Chapter  2  as 

?  n  ?  n  ? 

Kf  ,)       (  I      a„  cos  2tt  f  .t)  +  (  I     a.    sin  2tt  f  ±t )   ] 


(fi}  =  £  [    (  I      at  cos  27r  fit)  +  (   I   a^ 
1    n    t-1   t        x  t=l  t 


where  f.  =  i/n  is  the  frequency. 

2 
The  power  spectrum  p(f)  has  a  constant  value  2  a    over  the  frequency 

2 
domain  0=0.5  cycle  if  a  's  are  sampled  from  N(6,  °     ).   Consequently  the 

x*  f  a 

cumulative  spectrum  P(f)  =   w/     p(g)  dg  is  a  straight  line  if  plotted  against 

0 

f  (0  <  f  ^   0.5),  or 
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2 
P(f)/  a   is  a  straight  line  in  0  <  f  <  0.5. 

cL  - 

2  v 

P(f)/  a   can  be  estimated  by  2   I(fjj) 

a  i=l    J 

C(f.)  -   -1  S 

ns 

2  2 

where  s   is  an  estimated  variance  of  a  's,  i.e.   a 

t  a 

We  shall  refer  to  C(f .)  as  the  normalized  cumulative  periodogram.   The 
reason  why  the  cumulative  spectrum  is  used  rather  than  the  spectrum  itself 
is  that  the  Kolmogorov-Smirnov  testing  is  available  to  give  limits  with  the 
confidence  levels  around  straight  lines. 


Kolmogorov-Smirnov  testing  [9]  can  be  described  as  follows.   Let  F(x)  denote  the 
theoretical  cumulative  distribution  function  of  a  random  variable  (continuous) 
x,  and  F  (x)  an  empirical  cumulative  distribution  obtained  from  q  samples. 

D  =  max  |F  (x)  -  F  (x) |  takes  values  greater  than  or  equal  to   - — 
q         q  0°  2  /q~ 

with  probability  K  (A)  =  1-2   J   (~1)     e   ~       when  n  ->  °°.   Therefore 

v=l 
with  a  given  confidence  level  1-e  ,  F  (x)  should  be  confined  within  a  strip 

F(x)  ±  JLM     . 

As  for  the  case  in  this  section  F(x)  =  x  (0  <  x  <  .5)   and  q  =  — r"~ 

(n  even)  or  q  =  — —   (n  odd)  (*.*  Frequency  domain  is  (0,  .5)).   Therefore 

the  limit  lines  are  drawn  at  distances   K(e)/  ; —  above  and  below  the  line 

•q 

F(x)  =  x.   Some  values  of  k(e)  are  given  in  Table  6.1. 


e 

.01 

.05 

.10 

.25 

K(£) 

1.63 

1.36 

1.22 

1.02 

TABLE  6.1, 
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6.3  Correction  of  Model  Parameters 

In  some  cases  the  order  of  an  identified  model  is  correct  but  the  para- 
meters are  not  fit-ted  well.   It  is  possible  to  find  closer  values  for  para- 
meters by  analyzing  a  's. 

Suppose  that  the  correct  model  is 

<KB)  VdZt  =  9(B)  at 
but  that  an  incorrect  model 
<|>0(B)  VdZt  =0Q(B)  bfc 
is  used.   Then  b  will  not  be  Gaussian  any  more. 

bt  =  0o"1(B)  VB)  ^t  =  0o"1(B)  *o(B)  <f)~1(B)  0(B)  at  * 

That  is,  b   itself  is  an  ARMA  process.   Therefore,  by  identifying  b  in  the 
methods  described  in  Chapters  3  and  4,  for  example,  <f>'(B)  V  d  =  0(B)  a 
we  can  obtain  a  new  improved  model 

♦  '(B)  «.0(B)  Vd+d'zt  =  e0(B)  GT(B)  a^ 
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APPENDIX  1 
SOURCE  PROGRAM 
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THIS  FILE  PRINTED  ONi  02/13/72   AT i  05I19|53. 


str.iN 

file 

file 

INjEgE 
rEaL 
REaL 
%%%%%%%%%%%%% 
pR 


CARD( 

LIN 

R         I 


A 

lOCEOll 
VAL 
INT 
REA 
REA 

H 
X 
% 
% 
X 
% 
X 
% 
% 
X 
% 
% 
X 


KlNO  = 
ECKIN 
»  J.K» 
SIGM 
RRAY 
%%%%% 
RE 

UE  N 
EGER 
L  M 
L  ARR 
%%%%*> 
EGIN 


O.TITL 
D*7»MA 
P.Q.N. 
A  »  V  A  R  » 
ATCR. 
%%%%%% 
USID(Z 
*CO»L» 
O.CD.L 
EAN. VA 
AY  Z 
%%%%%% 


Eo6"N|jMBERS.")J 

XRECSlZE«22)| 

Ot  I  H  I 

MEAN.FPSILON, 

THETA.PHltOl 

%%%%%%%%%%%%< 

.N.MLaN.  ATCR, 

01 

,NJ 

RJ 

.Wt  U.ATCRCO 


.dlt.gama.zetm 

20]  > 
;xxxxxxxxxxxxx%xxx*xxxxxxxxxxxxxxxxxxxxt 

:,VAR»CD#L»D»wM 


1  » 


FnRMA 

Fjl(X 
F17(X 

ri?< 

Fl«(/ 

Fl6(X 
Fl5 


X 

X   R 

S    K 

%      L 
7 

?! 
X 
X 

% 

T  OUT 

lz»"o 

12»HN 
X12," 
/X5." 
"PART 
12»"V 
(X?.I 


PARA 
Z  «TIM 
W     ITIM 

N  INUM 
MEAN  lAVE 
VAR  IVAR 
ATCR  JAUT 
CD  IMAX 
0  «ORD 
L     IMAX 

THIS  PROC 
EGRESSIVE  A 

ESPECTIvELY 
IKELIHOOD  E 


METERS 
E  SERIES  TO 
E  SERIES  (DI 
RER  OF  DATA 
RA(,E  OF  Z   ( 
IANCE  OF  Z  ( 
nCORRFLATlON 
IMUM  LAG  OF 
ER  OF  HIGHER 
IMUM  LAG  OF 

EDURE  IS  FOR 
NO  MOVING-AV 
.THOSE  I N  T I  A 
STIMATES. 


BE  ANALYZED  (ORIGINAL) 

FFERENCErn 

POINTS (I.E.  Z) 

COMPUTED  IN  T"IS  PROCEDURE) 

COMPUTED  IN  T«IS  PROCEDURE) 

FUNCTION  (CnMPUTLo  HERE) 
ATCR 

DIFFERENCE 
PACT 

THE  INITIAL  ESTIMATES  OF  AUTO- 
ERAGE  MODELS,  THE  ORDER  IS  P  OR  Q. 
L  ESTIMATES  A*E  IjSED  IN  THE  MAXIMUM- 


pRnc 

VA 
IN 
RE 


EDURE 

LUE 

TEGER 

AL  AR 

%%%%% 

BEG 


RDER 

UMRER 

MEAN 

ORDER 

IAL  A 

ARIAN 

2»X9» 

REA 

DOU 

REA 

DOU 

TNT 

cu 

N»W1 
N» 

RAy 

%%%%% 

IN 
FDR 
F51 


OF  THE 
OF  DA 
VALUE 
"»X3»" 
UTflCOR 
CE  =  " 
Rll.4. 
L  SD 
BLE 

L  ARRA 
BLE  AR 
EGER 
%%%%%% 
RRELOP 

W» 
ATCRC 


DIF 
TAs" 
a  ". 
AUTO 
RELA 
•  XI. 
X  1  4  • 
J 

SUM, 
Y 

RAY 
I.J. 
%%f>% 
RINT 


FERENCE 
.X2. 14). 
XI  .  Rll. 
CORRELAT 
TIONS"// 

wil.a)» 
«ii.a)i 

SMI 

SQ. 
PACFCUL 
K.M. NN» I 
%%%%%%%% 

( ATCR.N. 


«  ".12/). 

4)» 
T0NSH.X3. 

)• 


PACCIOlLlI 

.1  IL1J 

NOX.  JNDX.QJ 

W)l 


01  J 

%%&%il%%%%%%%%%%%%%%%%%%%*%%%it%%%%%%%%%%%%%%l%%%%% 


MAT  OUT 

(X20."AUT0  CORRELATION  FUNcT I  ON"//  ) » 
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F52(X20»"V"»10(9("."),"V")). 

F53(X20»"A"#iO(9(","),"A")), 

F5«<  X  3 »  R 1 1 ,  4  ,  X? ,  1 2#  X  1 , "  I"  .  X50 .  "O" .  * (  "X"  )  )  . 

F55  <X3»RH.4.X2.I2.Xi,"I"»X*.*<"X"),"U«), 

F56(X3»R11.4fX2,l2»Xl."I".x50."0")» 

r57(X18."-l,0"(X7,"-.8"»X7»"-.6",x'»"-.4",x7,*.t2Mt 

X8."0",X8,nt2.t,X8,"t4".X8»  ". 6"»X8. ". 8". X?» "1 , 0" 

J  t 
F58(X20»"PARTIAL  AUTOCORRELATION  FACTION"//), 
F59(X10»"Q"»X5. "ERROR  BOUND"/). 
FA0CX9.I2.X5.Rll. 4)1 
INTEGER    I.J.KI 
REAL    TEMPI 

WRITE(LINE[SKIP  U)l 
IF  W«0  THEN 

WRITE(LINE#F51) 
ELSE 

wHiTErarNE:.F58)r 

WRITE(LINE[SPACE  2])l 
WRITE(LINE»F57)I 
WRITECLINE.F52)! 
FOR  I « »0  STEP  1  UNTIL  N  00 
BEGIN 

TEMPlsATCRri]! 
KI»50*TEMP| 
IF  K>0  THEN 

WRITE(LINE#F54.TEMP,I.K) 

ELSE 

IF  K<0  THEN 

WRITE(LINE.F55.TEMP.I.50+K.-K) 

Ei  r  r- 

WRITE(LINE.F56. TEMP. 1)1 
ENDJ 
WRITE(LINE.F53)| 
WRITE(LINE,F57)> 
WRlTE(LlNEtSPAcE   2I)j 
IF  Wao  THEN 
BEGIN 

WRITE(LINE.<X10."IF  THE  PROCESS  IS  MACQ". 
n).THEN  THE  AUTOCORRELATIONS  HIGHER  " 
"THAN   0     ARE  WITHIN  T«E  RANGE  OFl"//> 
)1 
WRITE(LINE.F59)> 
FOR  11=1  STEP  1  UNTIL  L  00 

WRITE(LINE»F60.I.SQtIJ)l 
END    ELSE 

BEGIN 

WRITE(LINE.<"IF  THE  PROCESS  IS  AR^t.thE  ". 

"PARTIAL  AUTOCORRELATION^  HIGHER  THAN", 
"  P  ARE  WIThIN  THE  RANGE  0Fl",X3.Rll,4//> 
•  SD)| 

ENO| 

ENOJ  %%      CORRELOPRINT  %% 

comment     this  procedure   is  to  calculate   THE  MEAN   v/ALUE. 

AUTncURRELATIONS  AND  PARTIAL  A.C.  PAC  IS  COMPUTED  IN 
OOU^LE  PRECISION  TO  GET  ENOUGH  ACCURACY,  I 
WRITE(LINE[SKIP  1))| 
FOR  I t ■ 1  STEP  1  UNTIL  N  TO 
W(!)i"'t!)l 
comment       higher  oifferencei 
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comment 


comment 


comment 


comment 


comment 


COMMENT 


Comment 
comment 


END 
AVERA 
SUM»« 
FOR  I 
SUM 
MEAN! 
■COVAR 

for  K 

BEG 


1  UNTIL  NN 
Wtlll 

NN  no 
1  UNTIL  CO  00 


STEP    1    UNTIL    NN-K    On 
:SUM*(W[ I3-MEAN)*(W[ i+K]"MEAN)| 
0    THEN 
>SUM/VAR  ELSE 


FUNCTIONi 


NN»«N| 

FOH  J I « 1  STEP  1  UNTIL   0  00 
BEGIN 

NNI«NN-1 I 

FOR  11=1  STEP  i  UNTIL  NN   00 
BEMN 

W[  1 3  |  eWC  1  +  1 1 
ENOIXSERIES 
I  XIORDER 
GEI 
0» 

l»l  STEP  1  UNTIL 
«=SUM+W[ I ]) 
aSUH/NNJ 
IANCEI 
1=0  STEP 
IN 

SUMlaO; 
FOR  I«»l 

SUMl 
IF  K  NEQ 
ATCR[K]  I 
BEGIN 

VAR»=SUMl 

ATCRtO) tail 
END! 
END» 

var»=var/nn; 

partial  auto  correlation 

PACCt 1  1  l"ATCR[ 1 II 

PACFClf 1 ] I'bnURLECATCRCll ) 
FOR  Jla2  STEP  1  UNTIL  L  00 
BEGIN 

PHKLiL)! 

SUMlaSMJ=OI 
FOR  M»  =  l  STEP  1  UNTIL 
BEGIN 

SUMlaSUM*PAcF[ J-l.M]  MUX  ATCK(J-M1» 
SM|a$M+PACF[ J-l.M]  MUX  AlCR[Mij 
ENOI 
PACF(J» J]  »»(OOUBLE(ATCRCJ])-SuM)/ll-SM)» 
PACCt Jl »=REAL(PACF[ J. J])l 
PHKL.K)  I 

FOR  K«=l  STEP  1  UNTIL  J-l  DO 

PACF[JiK]l»PACF[J-ltK]-pACF[J.J]* 
PACFt J-l.J-K]! 
ENOl 
TO  TEST.  THE  ASSUMPTION  THAT  THE 
SO»=«Ol 

FOR  Kj=1  STEP  1  UNTIL  L  00 
BEGIN 

SDI=SD+ATCR[K]**?I 
SQCK]  l>SQRT((l*Sn*SD)/NNl J 
END» 
THE  STANOARD  DEVIATION  SD  IS  C  ALC  llLATEU » 
TO  TEST  THE  ASSUMPTION  THAT  THE  PRUCESS  IS  AR(P), 
THE  STANOARD  DEVIATION  SO  IS  CALCULATED 
SD««1/(SQRT(NN) ) I 
WRITE(LINECSKIP  1])» 
WRlTE(LINE.</XiO.  "MESSAGES  FROM  US*OV/>)| 


J-l  DO 


PRUCESb  IS  MA(Q)I 
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WRITE(LINE»F17»N)I 
WRITE(UNE»F11.D>| 
WRITE(LINE»F12.HEAN)I 
WRlTE(LlNE»F16.VAR)f 
WRITE(LINE.F14)J 

WRlTE(LlNE»F15,0»ATCRtO)»PACCtO])| 
FOR  K««l  STEP  1  UNTIL  L  DO 

WRITE(LINE.F15,K.ATCR[K],PACF[K.K1  M 
CORRELOPRINT( ATCR»C0#O)  J 
C0RREL0PRINT(PACC#L»1)I 
WRITE(LINECSKIP  13)J 
END   PROCEDUREJ 
«XXXXXXSX*XXXXXXZXXXXXXX%XXXXXX£XX*XXXXXXXXX*XXXXXXXXXftXXXXXXXXXXXXXXXXXXXXXXXX? 
PROCEDURE  GAUSSLU(N»A»B»X) J 
VALUE      N! 
INTEGER    NJ 
REAL  ARKAY  Atl.l l.B.Xtl]) 

BEGIN 
FORMAT  OUT 

ri<X10#"MATRIX  WITH  ZERO  ROW  IN  DE^OMPUsE"/ ) . 
F?(X10»"SINGULAR  MATRIX  IN  DECOMPOSE."/)) 
INTEGER  ARRAY      PS(llN]l 
REAL   ARRAY        LUtllN.llNII 
FORMAT  OUT   FFl(XlOtRll«4«X10*Rll*4tX10fR11.4)f 
PROCEDURE  SINGULAR(WHY); 
VALUE      WHYl 


INTEGER    WHY! 
BEGIN 

IF  WHY»0  THEN 

WRITE(LINE»F1)J 
IF  WHYal  THEN 

WRITE(LINE»F2)I 
END!  %%    SINGULAR  XX 
PROCEDURE  DECOMPOSED)!   VALUE   NIINTEGER  N| 

COMMENT    COMPUTES  TRIANGULAR  MATRICES  I-  ANU  u  AND 
PERMUTATION  MATRIX  P  SO  THAT  LU,pA. STORES 
L-I  AND  U  IN  LU. ARRAY  PS  STORES  PtRMUTEn 
ROW  INDICES! 
BEGIN 

REAL  ARRAY  SCALEStl IN] I 
INTEGER    I»J»K»PVI! 
REAL       NR.PV»SIZE.MAX,MULTI 
LABEL    ELI 

COMMENT    INITIALIZE  PS»LU  AND  SCALES* 
TOR  Ilal  STEP  1  UNTIL  N  DO 
BEGIN 

pstni»n 

NRI'OI 

FOR  JlM  STEP  1  UNTIL  N  00 
BEGIN 

LUC  I . J  1  I  -  A  C  T • J  3  » 
IF  NR<ABS(LlJtI»J])   THEN 
NR|aABS(LUtl»J])| 
END  I 
IF  N<*  NEQ  0  THEN 
SCALESCI] |«1/NR  ELSE 
BEGIN 

SCALFSt II laOt 
SINGULARC0)» 
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END1 
END  I 

COMMENT    G 

P 

FOR  KJ«1  ST 

BEGIN 

MAXIMO 

FOR  H 

BEG 

S 

IF  MAX 


IF 


IF 


END 

MAX 

BEG 

S 

G 

END 

PVI 

BEG 

J 

P 

P 

END 

PVi*LU 

FOR  II 

BEG 

L 

I 


EL« 


ENDI 
IF  L 


END 

UtPS 

SIN 

OSE 


AljSSIAN  ELIMINATION  wiTH 
ARITIAL  PIVOTING! 
EP  1  UNTIL  N-l  DO 

J 

■K  STEP  i  UNTIL  N  DO 
IN 

IZEI»ABS(LU(PS[I1.K)  )*SC  ALES!  PSU  H  > 
<SIZE  THEN 
BEGIN 

MAXI.SIZEI 

PVI  lal 

ENDI 
I 

■0  THEN 
IN 

INGULAR(I )l 
0  TO  ELJ 
i 

NEQ  K  THEN 
IN 

l=PS[K]l 
S[K]  I"PS[PVI]I 
SCPVI JI"JI 
I 

t  PSt  K  3  »K  i  I 

=K+1  STEP  1  UNTIL  N  D" 
IN 

urpsr  1 3  »k3  i  =  mijlt i*Lnr^St  I J  »K3/pVl 

F  MULT  NEQ  0  THEN 

FOR  JlsK+1  STEP  1  UNTIL  N  DO 

LU(PS[I]iJ]l«LU[P^(I]».l]-MULT* 
LUtpi>CK]»j]| 
) 

CN]  tNl»0    THEN 
G|jLAR(l)| 

%% 


endi   %%  decompi 
procedure  .solve   j 

COMMENT   SOLVES  AXoB  USING  LU  FROM  DECOMPOSt-t 
BEGIN 

INTEGER 
REAL 

FOR 
BEi 


I.J* 

DOT! 
I«»l 

:gin 

DOT 
FUR 


STEP  1  UNTIL  N  nO 


ENDI  *% 


Xtl 

ENDI 
FOR  II»N 

BEGIN 
OUT 
FOR 

XII 

ENDI 
SOLVE  %% 


IbOI 

Jlai  STEP  I  UNTIL  1-1  DO 
DOTl«nnT*LU[PStIJ»JJ*X[J]l 
]|=B[PStI]]-DOTI 

STEP  -1  UNTIL  1  DU 

1=0  J 
J  t ■  1*1  STEP  1  UNTIL  N  DO 
00Tl=0nT  +  Ll)[PSCI3.ji*X[J]> 
] i»(X[  I3-00T)/LUtPS[I  J»I1» 
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xxxx 

X 
X 


MAIN  XXX 


ENDI 

tgsszxxxxxxxxxxxxxx 

proceoure  us 
value  n. 
integer 
real  mea 

REAL  ARRA 
X  INPUT  P 


N 

MEAN 

VAR 

ATCR 

P 

PHI 

0 

THETA 

SIGMA 

EPSIL 


XXXXXXXXXxXxXXXXXXX 

REGIN 
rORMAT  OUT 
FfKXIO," 

RU.4 
FF2(XlO." 

X2.I2 
Ff3(X10." 
FF33(X10. 
Ff^CXIO." 
Ff5(X10." 
REAL 

INTE 
REAL 
LABE 
XAR- 
X 

FOR 
RE 


EN 
X 

GAUS 
PHI  t 

FOR 


IF  N 


ELSE 

IF  N 

BE 


END 
XX  G 

%%%%% 
PE(N. 

P.  Q» 
N,P*Q 
N.VAR 
Y  ATC 
ARAME 
|NUM 

•  AVE 
IVAR 
IAUT 
lASS 
IPAR 
lASS 
IPAR 

•  w  A  R 
OMlCR 
%%%%% 


\  THEN 
IF  All 

X 
ELSE 

S 

NEQ  0 

GIN 
DECOMP 
SOLVEI 

I 

AUSS-LU 

xxxxxxx 

P.Q» 

MEAN* 
I 

.SIGMA* 
R. THETA 
TERS 
BER  OF 
RAGE  OF 
IANCE  0 
OCORREL 
UMED  OR 
AMETERS 
UMEO  OR 
AMETERS 
IANCE  0 
ITERIDN 
%%%%%%% 


#11  NEG  0  THEN 
ti)l"B[l]/Atl.lJ 

INGULAR(I) 

THEN 

OSE(N)l 

XX 

%%%%%%%%%%%%%%%%% 
mean.var.atcr.  p 

VAR.ErstLONI 

EPSILON) 
■PHI |*ll 

w 
w 

F  W 

ATIONFUNCTION  OF 
DER  OF  AUTn-REGRE 

FOR  AUTO-REGRESS 
OER  OF  MOVING-AVE 

FOR  MOVING-AVERA 
F  GAUSSIAN  TERM  ( 

FOR  CONVERGENCE 


XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXZ'* 

HI»THElA»SlGMA#EoSlLON)J 


DATA»"»I3.X5."MEAN»".R1 
".I2.X10 


Nn.OF 

/). 

ORDER  OF  AR-PROcESS 
/). 

»EPSIL0N«".Rll,4/). 
AR-PROCESS  PARAMETERS". X15»"M 
PH 1 1 " » 1 2 .  "  1  " .  X 1 ,  R 1 1 ,  4 .  X20.  "TH 
ARRAY  AUIP.HP].  RtltP]. 
TtitQ+l.UQ  +  ll.TAU. 
GFR    I.J.K.L.Mf 

MAXl 
L  LOUT  I 
PARAMETER  ESTIMATED 


Iisl  STEP  1  UNTIL  P  00 

GIN 
R[I]I-ATCR[Q  +  U*VARI 
FOR   Jt « 1  STEP  1  UNTIL  P  DO 

AII.J1 l»ATCR[ARS(Q*I-J)]*VAR| 

Dl 

SLU(P»A.B.PHT)I 

I  |31  STEP  1  UNTIL  P  DO 


W 

SSIVE  MODEL 

IVE   MUDEL  fTO  BE  ESTIMATED) 

RAGE  MUDEL 

GE  MODtLCTU  RE  ESTIMATED) 

TO  SE  t-STlMATED) 

IN  NEwlON-KAPHSON  S  MFTHOD 

XXXXXX»XXX»XXXXXXXXXXXXXXXXXXXX* 


1.4.X5»"VANlANCE»"» 
."ORDEM  OF  MA-PROCESS". 


A-PROC^-SS  PARAMETERS"/). 
ETA[".A2»"l«»Xl»Rll.4)l 
PHTCllfJ. 
C.F.H[llQ+l]| 
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PHI  III l"PHTC III 

X 

%    MA-PROCESS   PARAMETER  ESTIMATE 
X 

X(l)  MODIFIED  COVARIANCE  SEQUENCE 
X 

FOR  Il»l  STEP  1  UNTIL  Q+l  DO 
RE6IN 

IE  P  EU  0  THEN 

CCI1  l«ATCRCI-l]*VAR 
ELSE 
BEGIN 

MAXlaOI 

FOR  Ji»0  STEP  1  UNTIL  P  DO 

FOR  Kt»0  STEP  1  UNTIL  P  UO 

MAXI«MAX*PHI(J]*PHIlK]*ATCRtABS(J*I 
•K-l)] J 
Ct  I  I  l»MAX*VARI 
ENDl 

END!   X  I  LOOP 
TAUtl] »=SQRT(Ct 1])» 
FOR  Iis2  STEP  1  UNTIL  0*1  nO 

TAUt I  1 »=0J 
XNEWTON-RAPHSON  METHOD 
FOR  Lisl  STEP  1  UNTIL  10  On 
REGTN 

EOR  1**1  STEP  1  UNTIL  Q+l  DO 
BEGIN 

FOR  J I » 1  STEP  1  UNTIL  0*1  DO 
BEGIN 

IF   I+J  LEQ  Q+l  THEN 
BEGIN 

IF  J  GEO  I  THEN 

T[I#J3lsTAU[I*J-n  +  TAUtJ-I*l3 
ELSE 

T[ItJllaTAU[I*J-lJJ 
ENn 
ELSE 
IF   J  GEO  I  THEN 

T(I»J] |«TAU[ J-I+ll 
ELSE 

T  C I  •  J 1 1  ■  0 1 
ENOI  X  J  LOOP 
MAX»=0! 
FOR  JIM  STEP  1  UNTIL  0-1*2  Ou 

MAXJaMAX  *TAUt Jl*TAUtI*j"l Jl 
FCI1I«MaX-C[ III  ' 
END»   XXXI  LnOP 
GAUSSLU(Q+1»T.F»H)I 
FOR  Il»l  STEP  1  UNTIL  0*1  DO 

TAUt  II  «  =  TAUt  I  1-Ht  ID 
MAX»=ABS(Ftl J)| 
FOR  J I =2  STEP  1  UNTIL  Q+l  DO 
IF  MAX<ABS(F[ Jl)   THEN 
MAX:=ABS(F[ J])l 
IF  MAX  LEO  EPSILON   THEN 
GO  TO   LOUT! 
ENDi   X  LOOP  L 


LnUTl 


FOR  Ii*l  STEP  1  UNTIL   Q  DO 


-57- 


THETAtlJlB-TAUn-Ml/TAUCHl 

IT  P>0  THEN 
BEGIN 

maxi«oi 

FOH  I l«l  STEP  1  UNTIL  P  00 
MAX|BMAX*PHIt I] I 
THETA[03t«MEAN*(l-MAX)| 
END    ELSE 

THETAC03i»MEAN| 
IF  0>o  THEN 
BE<MN 

SlGMA|«TAUtl  J*TAll[13l 
ENO 
ELSE 

beGin 
maxi*oj 

FOR  Il*l  STEP  1  UNTIL  P  DO 

MAXlaHAX*PHI[n*ATCRtIll 
SIGMA|bVAR*(1-MAX)| 
ENOl 
WRITE(LINE[SKIP  1])| 

WRITEfLlNE.<X10#wMESSAGES  FROM  USPE"//>JJ 
WRITE (LlNEtFFl.N.MEAN.VAR)i 
WRITE(LINE.FF2»P#Q)» 
WRITEfLlNE.FF3.SIGMA)J 
WRITE(LINE.FF33.EPSIL0N)I 
WRITE(LINE.FF4)1 
IF  P>Q  THEN 
BEGIN 

FOR  I»*0  STEP  1  UNTIL  Q  00 

WRlTE<LlKiE.FF5»I.PHICI3»I.THET*[l3)| 
FOR  Il"Q*l    STEP  1  UNTIL  P  00 

WRITE(LINE»FF5»I,PHItI3»I»0)| 
ENO    ELSE 
BEGIN 

FOR   1 1=0  STEP  1  UNTIL  P  00 

WRITE<LINE»FF5»I,PHIU3#I.THEIA£IJ)I 
FOR   II«P+1    STEP  1  UNTIL   Q  DO 

WRITECLINE»FF5.I,0.  I.THETAtHM 
ENOl 
ENO>  %    USPE 
%%%%%**%**%*%%%%%%%%%%%%%%%%%%%%%%%%%%%%l%%%%%%%%%%%%%*t%%*%%%%%%%%%%%%%%%%%%%%i 
PROCEDURE    RFSlDUAL(W.N»PHI#P.THETA#Q.MEANtA)J 
VALUE  N»PtQ*MEAN; 
INTEGER     N.P.Qj 
REAL    MEAN! 
REAL  ARRAY   PHI # THET A . w , A [ *  ]  | 

%%%%%*%%%%*%%%%%%%%%<K%%%i%%%%%%%X%%%%%%%%%%%%%%%%%%%%%*%%%%t%%%%%%%%ltt%%%%%%%%} 

BEGIN    COMMENT 

ThIS  SUBROUTINE  COMPUTES  THE  RESI0UAL  TErmS ( GAUSSI AN 
TERMS)  UF  TIME  SERIES  BY  BACK-FURECAS • I NG  AND  FORWARD 
FORECASTING. THIS  SUBROUTINE  IS  TO  TE  CALLED  IN  MLES(MAy 
IMUM  LIKELIHOOD  ESTIMATES)  AND  USNLCN^N  LINEAR  APPROxI 
MATI0N)> 
INTEGER  I.J.K; 
REAL     TMP.TMQI 

REAL  ARRAY   E  .  X t -1 0-M AX( P . 0 ) l N  +  Q-P ] I 
FOR  Il«-10  STEP  1  UNTIL  N  DO 
BEGIN 

At  I  )  1  =  01 
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IF  I  GEQ  1  ANO  I  LEQ  N  THEN 
BEGIN 

X[ U «»w{  I3-MEAN) 
END 
ELSE 

Xt IJ iaO> 
ECHiaOi 
ENOl 
RACK-FORECASTING  FOR  F'S  AND  *'S. 


TOR  I  I-N-P  STEP  -1  UNTIL  1  DO 
BEGIN 

TMPnTMQlaOl 

FOR  J»=l  STEP  1  UNTIL  P  00 

TMPJaTMP  +  PHlt J  1 *X[ J*I]  I 
FOR  Jl  =  l  STEP  1  UNTIL  0  nO 

TMQI»TMO+THETACJ]*ECJ*I]» 
EIDUMI  1-TM-P*TMQ) 
ENO) 
TOR  1 1 =0  STEP  -1  UNTIL  -10  DO 
BEGIN 

TMPt*TMQteOJ 

FOR  Ji«l  STEP  1  UNTIL  P  DO 

TMPl«TMP*PHlt J1*X[ J+IJI  . 
FOR  Jl«l  STEP  1  UNTIL  Q  DO 

TMQl.TMQ+THFTA[JJ*Et J*I]' 
Xtll «=E[ Il+TMP-TMQI 
ENOJ 
%  FORWARD  FORECASTING  FOR  A'S 

FOR  Ki«-10  STEP  1  UNTIL  N  00 
BEGIN 

TMPlsTMQlxOI 

FOR  J t  =  1  STEP  1  UNTIL  P  DO 

TMPlaTMp+PHlC J]*X[K-J] I 
FOR  Ji=l  STEP  1  UNTIL  0  DO 

TMQt»TMQ*THFTA[J3*A[-J*KJJ 
ACK]  l*XtKl-TMP*TMQj 
END! 
FNn  OF  RESInilALl 
%%%%i%%it%%%%%%%%%%%%%t%%%%%%%i%t%%%%%%%%%%%%%%%t%%%%^%%%^%%%%%%%%l%t%%%%%%%%%i 

real  procedure   mles(w.n.phl»p.thfta»q.meani ls'» 
value    n»p»Q»mean» 
real  mf.an.ls> 

INTEGER    H,P'Ql 
RfAL  ARRAY   PHI.THETA.WC*]  | 

%%%%t%%%%%%%%%%%%%%%%%ii%%%%i%%%%%%%%%%%%%%%%%%%%%%%%*%%%*%%%%%%%%%%%%%%%%%%%%y 
REGIN 

COMMENT     LEAST  SQUARE  FUNCTION  IN  CUMPUTFD  FOR 
THE   RIVEN  RESIDUAL  TERMSJ 
REAL  ARRAY   A  [  -  1  O-UA  X  (  P  #  Q  )  !  M  ]  I 
REAL      TMp; 
INTEGER    I; 
RESIDIIALCW,N»PHI»P.THETA.Q,MEAN,A)J 
ILEAST-SQUARE  FUNCTION 

TMPj=n» 

FOR    Iia-10    STEP    1    UNTIL    N    r>0 

TMP«*TMP  +  A[ I]*At  III 
LSlaTMPJ 
'    ENDJ  SSMLES 

%%%%i%%%i%r%%f%%%%%%%%zt%%%i%%i%%%%%%%%z%i%%%%%%%%i%%*%%%*%%%%%%%%%%%%%%%%%%t%^. 
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pRflCrDURE       DGCKCA.N.PQ.T1.ITER)! 
VALUE         N.PQ»T1.ITERJ 

INTEGER   n.hq.iteri 

real   ti» 

real  array  A[*]i 

t%%%%*%%tX%%%%%%l%%t%%%%%%%%%%%%%tt%%%t%%%t%%%%%%%%%%l*%%%%%%tl%%%%%%t%%%%%%%%%1 

REGIN 
REAL      MEAN.TMPfTOtCHISQl 
INTEGER     I,J.K> 
TMPiaOl 
FOR  I|a-10  STEP  i  UNTIL  N  DO 

TMPl*TMP+ACI3l 
MEANIaTMP/NI 
CHISOl»OJ 

FOR  Ki=0  STEP  1  UNTIL  ?0  00 
BEGIN 

TMPt«o» 

FOR  lJs-10  STEP  1  UNTIL  N-K  On 

TMPl«TMP+(ACI]-MEAN)*(AtI*K)-HEAN>l 
NEQ  0  THEN 
CHISQ«aCHIS0+(TMp/T0)**2 


IF  K 
ELSE 


ITER.T1* 


TO«»TMPl 
END  OF  K  LOOP* 
CHISQ|*CHISQ*NI 

WRITE(LlNE.<X«tI3.X59»Rll.4tX2.Rll.4»X7»I3  > 
CHISQ.K-PQ-D1 
ENO  OF  DGCK1 

PROCEDURE    MlPRINT(Z»N»PHI.P#THETa»Q»MEAN.OLT)' 

value        n.p»q»mean.dltj 
Integer  n.p.qi 

real        mean»dltj 
real     array        ZpPhi.thetac*]  i 
%%%%%**%%%  Jt%%%%%%%%%%%%i%%%%%%%%%%%%%%%%%%%%%%i%%%%%%%*%%%t%%%%%%%%%t%%%%%%%%%%t 
% 


N 

Z 

PHI 

P 

THETA 

Q 

MEAN 
OLT 


% 

% 
% 

% 
% 

% 
% 
% 

% 
% 
% 

% 
% 
% 

% 
% 
BEGIN 

FORMAT 


NUMBER  OF  DATA 
TIME  SERIES 

AUTO-REGRESSIVE  PARAMETERS 
THE  ORDER  OF  AR  PARAMETERS 

MOVING  AVERAGING  PARAMETERS 

THE  ORDER  OF  MA  PARAMETERS 

THE  MEAN  VALUE  OF  Z 

SMALL  CHANGES  OF  PARAMETERS 


THIS  PROCEDURE  PRINTS  OUT  THE  SUM  OF  T HE  SQUARES  FUNCTION 
WHEN  THERE  ARE  1  OR  2  PARAMETERS. 

PARAMETER 
11*11  VALUES 


DTMt 
DIM. 


21  VALUES  AROUND  THE  PRE-F^TI MATED 

AROUND  THE  PRt-ES T I  MATED  PARAMETERS. 


OUT    FM1(H(X2.R9.  ?)//). 
FM2C11CX2.R9. ?)//). 

FM4(R'?,2."*".ll(X?.R9.2)/.A9."*V,X9,"*"). 
FM3(X10.11(X2.p9,2))l 

CTHETA[0|Q).cPHlt0lP].LS["5l5J| 


REAL  ARRAY 
REAL    A.b#PA) 
INTEGER    I»J.KI 
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WRITE(IINE[SKIP  11)1 

WRlTE(LlNF.<XlO. "MESSAGES  FROM  MLP«INT   (THE  SUM" 
•  "  OF  THE  SQUARE  FUNCTION"/'/:*)* 
IF  P*o«l  THEN 
B  E  G  T  N 

FOR  Jl=-5  STEP  10  UNTIL  5  00 

regin 

for  i«»"5  step  1  until  5  do 

BEGIN 

IF    Psl     THEN 

CPHI[ 1] js    PHI[ 11*0LT*( UJ) 
ELSE 

CTHETACl]|s    THETAtl]*OLT*(I*J)l 
MLES(Z»N»CPHIfP»CTHETA#QtMEAN»LSCl])l 
end; 

IF    P=0    THEN 

PAl*THETA[ 1 ) 

ELSE 

PAlsPHIt 1 ] » 
WRITE(LINE»FM1.F0R    Il«-5    STEP    1    UN.l.IL    *    00 

PA*(I+J)*DLT)» 
WRITE(LINE»FM2.F0R    II.-5    STEP    1    UNUL    S    00 

Lstini 

WRlTE(LlNEtSPACE  23)1 
EN01 
END   ELSE 
REGIN 

IF  PaO  THEN 
BEGIN 

AlsTHFTAC  1  3  I 
BI*THETA[23  J 
WRITE(LINE.</X20,"A«THETA[ll"»X5»"RxTHETAC2]" 

END   ELSE 
IF  PM  THEN 
BEGIN 

A  I  =PH I C 1 3  I 
BI=THETAt 13 » 

WRITE(LINE»</X20."A*PHltl3,,.X'»"BsTHETA[l)"//> 
) 
ENO  ELSE 
BEGIN 

A-laPHIC  111 
B«aPHl[23» 

WRlTE(LINE.</X20,"A*PHItl3,,.X''»"BePHlt2]"//>)l 
ENDS 
WRITE(LINE.<X115."A"/>)J 
WRITE(LINE»FM3.F0R    I»a-5    STEP    1    UNUL    * 

00    A*I*DLT)J 
WRITE(LINE»</X9,121{"*")>)| 
FOR    Iia-5    STEP    1    UNTIL    5    00 
BEGIN 

IF    P=0    THEN 

CTHETA[23»=  THET A [ 2 l *DLT* I 
ELSE  IF  Psl  THEN 

CTHETAt 13  »sTHETA[13  +  I*DL' 
ELSE 

CPHI[2]|«  PhIC21+DLT*II 

FOR  J I  =  - 5  STEP  1  UNTIL  5  DO 
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procedur 
value 

REAL 

INTEG 

REAL 

% 
X 

% 
% 
% 

% 
% 
% 
% 

% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
XT 

begin 
Forma 
Fr^cx 

FF6<" 


Ff5(x 

rea 

INT 
REAL 

LAR 
REA 


BEGIN 

IF  P"0  THEN 

CTHETA[1JI»  THtTAtll*OLT*J 
ELSE 

CPHItlJl«PHI[l J*J*UlTI 
MLES(Z#N.CPHl»P»CTHtTA»U.MEAN.LS[J]) 

I 

ENni 

WRITF(LINE»FM4.R*I*0LT»FUR  JI.-5  STEP  1 
UNTIL  5  00   LS[JJ)| 
ENDI 
WRITECL!NE.</X3."B">)I 
ENDl 

WRITE(LINEtSKIP  U)l 
ENDl  %%    SUM  OF  THE  SQUARES  PRINT  %% 

%%%t%%*txz%%%%%%%%%%%%%%%%%%%%t%%%%%%%%*%t%*%%%%xt%%%xt%%%%%%%%%? 

E   MLQT(WtN.PHIiP.THETA»Q»MEAN.EPSILON.UElTA,iM#FACTOR, 
LAemDA.LMAX) , 
N.P,«»MEAN#EPSILnN»OELTA»IM.FACTOR#LAMfloA,LMAX| 
ME  AM, EPS  RON. DELTA. FACTOR .LAMBDA . LMAX* 
ER    N.P»Q.IMI 
ARRAY   PHI.THETA.Wt*]! 
%%%%%%%%%%%%%%%%%%%%%%%%%%X%%%X%%%%tX%%*%X%*%%%%%l%Xl%%X%l%%%%%% 


W    =   TIME  SERIES 

N       NUMBER  OF  DATA(W) 

PHI    AR  PARAMETERS 

P       URDER  OP  PHI 

THETA   MA  PARAMETERS 

Q       ORDER  OF  THETA 

MEAN    AVERAGE  OF  * 

EPSILON  MAXIMUM  ERROR  OF  CONVERGENCE  OF  PARAMETERS 

DLT     DELTA  USED  TO  CALCULATE  THE  DERIVATIVES 

IM     MAXIMUM  NUMBER  OF  ITERATION 

FACTOR  PARAMETER  TO  MODIFY  LAMRDA  DEPENDING  PN  THE  CONVERGENCE 

LAMBDA  PARAMETER  TO  MODIFY  THE  DIAGONAL  ELEMENTS  IN  MARQUARDT'S 

ALGORITHM 
LMAX   MAXIMUM  OFLAMflDA 


THIS  PROCEDURE  CALCULATES  THE  PARAMETERS  ITERATlVELY, 
CONVERGENCE  IS  BETTER  THAN  GAUSS-NEWTON ! S  MtTHODC USNL) . 
THE  PROGRAM  IS  ALMOST  IDENTICAL  TO  USNL  EXCEPT  THftT 
DIAGONAL  ELEMENTS  ARE  MODIflED  TO  GIVE  OETTtR  CONVERGENCE. 
HE  DERIVATIVES  ARE  CALCULATED  AT  FIRST. 

T         OUT 
105.2(.X2.R9,2)/>» 

*    OF    ITER.  ".X9."PHI"»X20. "THETA".  X18»"i>UM    Ur    SQUARES". 
X2. "CHI-SQUARE". X2.HQEG.    OF    FREE . " . X2. "L AMRU a« , 
X4."MAXBETA"//). 

12."PHlt"»I2»"]"»Xl.Rll.A.X04»"THETA[".J2."J"»Xl»Rll.<ni 
L       TMP.DELTAINV.MAC.LS.LSOI 
EGER  I.  J.K.PQ. ITrRl 

ARRAY    A»AlflO-Q|N  +  P*Q]»X[liP*Q«-10lN*»'*01»fllllP*0» 
llP+QJ.F»BETAtllP*Q]l 
EL  L2.LSICK.LDUTI 

L    ARRAY  D[llP*Q]J 
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mles( 

WRITE 
WRITE 
WRITE 
DELTA 

PQl«P 
RESID 
FOR 
REG 


W.N, PHI 
(LlNEtS 
(LlNE.< 
(LINE.F 
TNV I ■ 1 / 

UAKrt.N 
ITER  I al 
IN 

FOR  !•■ 

BEGIN 

IF 

EL 

RE 
IF 

EL 

F 

ENDJ 
FO 


.P.THETA#Q.MEAN#LSO)l 

KIP  11)1 

M5»"MESSAGES    FROM    MLQT"//>)I 

F<S>» 

DELTA) 

•PHI»P»THETA»Q,MEAN»A)I 
STEP  1  UNTIL   IM  DO 

1  STEP  1  UNTIL  PQ  DO 


I 

SE 

SI 

I 

SE 
OR 

% 

R 
BE 


LEQ  P  THEN 

PHICI  Jl«PHI[M  +  DELTA 

THETAtl-PJ  l«THETA{I-P]  +  Dk-LTA» 
DUAL(W.N,PHI.P.THETA,Q.MEAN»Al)| 
LEQ  P  THEN 
PHIlrl'aPWUr  r-UELTA 


THET 
Jl*- 
XCI, 
OF  C 

I«"l 

GIN 

FOR 
BE 


FO 
FO 


EN 
DC  I J 

TMPI 
FOR 

Fen 

D  OF 

Jl«I 
IF  I 


L2t 


FOR 


A[I-P]  »* THET At  I-P]-D*-LTA» 
10  STEP  1  UNTIL  N  DO 
JJts<A£j]-AlCJ))*DELlAlNV| 

OMPUTING  DERIVATIVES 
STEP  1  UNTIL  PQ  DO 

Jl *  I  STEP  1  UNTIL  PO  DO 
GIN 

TMPIbOI 

FOR  KIb-10  STEP  1  UNTIL  N  DO 

TMPl»TMP*XtNKi*X[  J.KJI 
BtI«J3:«TMPj 
IF  I  NEQ  J  THEN 

Br  j,n  ibtmpj 

Dl  %    OF  J  LOOP 
I «SQRT(BC  X • 1 3  > I 

=  01 

KJ«-10  STEP  1  UNTIL  N  00 

T"PiaTMP*Xt  IiK]*ACkJ  J- 
t  ■  T  M  P  /  0  [  I  ]  I 
I  LOOP! 

STEP  t  UNTIL  PQ  DO 
STEP  1  UNTIL  PQ  DO 

NEQ  J  THEN 

BtI»J]|"BtJ»I3l*BtI»Jj/tDtI3*D[J3)l 


DO 


I  1 ■ 1  STEP  1  UNTIL  PQ 

BClt  I)  lol+LAMRDAI 
GAUSSLU(PQ,B,F#BETA)I 
FOR  1 1 » 1  STEP  1  UNTIL  PQ  DO 

BETAt  I]i=BETAm/DtI]l 
FOR  H*l  STEP  1  UNTIL  PQ  DO 

IF  I  LEQ  P  THEN 

PHimi*PHim*BETAin 

ELSE 

THET At I-P] «= THET At »P]+RtTA[ 1] I 
MLES(W»N»PHI»P#THETA.Q#MEAN»Li>>l 
IF    LS    <    LSO    THEN 
BEGIN 

LSOicLS) 
RESI0UAL(W,N.PHI#P»THETA.Q,MEAN»A)I 
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LOiiTi 


TMPt.oi 

FOR    Ki«-iO    STEP    1    UNTIL    N    00 

BEGIN 

TMP|BTMP*ACK]*AtKll 

ENDl 

OGCK(A»N»P*Q»TMP»ITER)l 
IF    P>Q    THEN 

BEGIN 

FOR    ll«0    STEP    1    UNTIL    Q    00 

WRlTE(LlNE.Fr5»I.PHItIJ.I.THETA(IJ)» 
TOR    It«Q*l  STEP    1    UNTIL    P    00 

WRITE(LINE#FF5.I,PHICIJ»I.0)| 
ENO         ELSE 
BEGIN 

FOR      I « - 1    STEP    1    UNTIL    P    00 

WRITE(LINE.FF5.I.PHI[I].I.THElAcn>| 
FO.R      II*P*1         STEP    I    UNTIL      Q    0X1 

WRITEaiNE»FF5»I,0#I.THETA[Il>l 
ENDl 

MACl«ABS(8ETAtlJ)l 

FOR  11*2  STEP  1  UNTIL  PQ  00 

MACl»MAX(MACtABS(8ETACI] >>l 
IF  MAC<EPSILON  THEN  GO  TO  LOu'l 
WRITE(LINE»FF4.LAMBDAtMAO) 
LAMBOAI«LAMBOA/FACTOR) 
ENO    ELSE 
BEGIN 

FOR  1 1 » 1  STEP  1  UNTIL  PQ  00 
IF  I  LEQ  P  THEN 

PHX(X3l"'PHlCX]aBETA[  j] 
ELSE 

THETA[I]l«THETAtn-RETA(ni 
LAMBOA«aLAMROA*FACToRl 
IF  LAMBOA>LMAX  THEN  GO  T"  LSiCKl 
GO  TO  L2I 
ENOI 
ENOl 
LSICKI 

IF  P>0  THEN 
BEGIN 

TMPl-OI 

FOR  1 1  - 1  STEP  1  UNTIL  P  00 

TMP|«TMP*PHICI]| 
THETA[0J|«MEAN*(1-TMP)| 
END  ELSE 
THETACO]laMEANi 

WRITE(LINE.</X35.-THETA[0]-".R11.4///>»THETA[03)I 

WRITE(LINE»<X123."*">)I 

WRITE(LINE#<X123. «*">)! 

WRITE(LINE.<X10.WN0  MnRE  CONVERGENCE. CHFCK  IF  "• 

"SOLUTIONS  ARE  REASONABLE  OR  NUT  .  "»52<  "*•»>'>  )  I 
WRITE(LINE»<X5."IF  MAXBETA  IS  VERY  SMALL  OR  ThE  •», 
••CONVERGENCE  OF  PARAMETERS  IS  OBVlUUSL*  GOOD."* 
•♦yOU  ARE  GETTING  RIGHT  ANSWERS.   0N  TH&  OTHER  ". 
"HAND. IF  MAXBETA  IS  STILL  LARGE  AFlER  Max  " . 
MIERATION  OR  IF  THE  SUM  OF  SQUARE*  IS  NOT  ". 
••APPROACHING  THE  MINIMUM  VALUE.  YOU  SHUuLO  ". 
•.RE-EXAMINE  THE  DATA.  THIS  MESSAGE  COMts  WHEN  ". 
"LAM30A  REACHES  OR  EXCEEDS  "» I 4/>,LMAX> I 
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End  of  mlqt 
r^aoccard./ 

BpGIN 

REAL  ARRA 
LABEL 

FOR  lis 

READ(CA 
FINISH! 

WRITECL 
FOR  It* 

W 

USIDCZ, 
WRITE(L 
WRITECL 
FOR  It. 

W 


I 

%*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*%%%* t%%%%l%%%%%%%%%%%%%t  i 

%%%%%%%%%%%*%%%*1,t%%%%%%%%%%%Z%%Zl%%i 
MAIN  PROGRAM        %%%%%%%%%%%*l%%*%%t,%%%%%l%l%xli,i%i%1<i 

%l%%%%%%%%%%%%%%*%%*%%%*%%%l%%%l%Z%*%%tt%%l%i 

T.EPSILON.GAMA.ZETA" 


,N»P.D#0»  IM»t)L' 


Y  W.Ztl  lN*iO]pA[-15lN*10]l 

FINISH! 

0  STEP  1  UNTIL 

RU./.FOR  Jtal 


N  01 V  7  no 
STEP  I  UNTIL  7  DO  Zl7*I  +  J]  HFINI'SH]  J 


END. 


USPE(N- 

MLPRINT 
MLQT( N»N- 
ENn» 


INE.< 
0  STE 
RlTtC 

N»MEA 
INECS 
INE»< 
0  STE 
RITE( 

0»P»Q 

(  tit  N- 
D»PHI 


X10."lNPlj 
P  1  UNTIL 
LINE»<X2, 
DO    Z  C 1*7  + 

N» ATCR# VA 
KIP  1])J 
X10."M0DI 

P  1  UNTIL 
LINE.<X2, 
DO  W[I*7* 

.MEAN#VAR 

D»PHl»PtT 
.P.THETA, 


T  DATA  (TIME  SER I ES *"///>  )  I 

N  DIV  7  DO 
7(F8.5.X3)>»F0R  Jlal  STEP  1  UNTIL  7 
JDI 
R»?0,2o»D»W)| 

FIED  TIME  SERIES"///>)J 

N  OIV  7  00 
7(r8,5.X3)>.F0R  Jlal  STEP  1  UNTIL  7 
J3)> 

.ATCR.   PHI,THETA#SAGMA»EPSIL0N)I 
HETA»Q»MEAN#OLT)l 
Q»MEAN,EPSIL0N»DLT»*0.2.0..01»100)| 
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APPENDIX  2 


ARMA   (1,    0,   1)   Model 
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Input  data  (time  series) 


17,00000 

16 

.60000 

16 

.30000 

16, 10000 

17,10000 

[6.90000 

16, 

80000 

lf.toooo 

17 

.10000 

17 

,00000 

16,70000 

17.40000    1 

['.20000 

17, 

40000 

17,40000 

17 

.00000 

17 

,30000 

17,20000 

17.40000    1 

1&.80000 

17, 

10000 

17,40000 

17 

.40000 

17 

,50000 

17.40000 

17,60000    1 

['.40000 

17, 

30000 

17.00000 

17 

.80000 

17, 

.50000 

18. 10000 

17,50000    1 

', 40000 

17, 

40000 

17.10000 

17 

.60000 

17, 

,70000 

17,40000 

17.80000    1 

['.60000 

17, 

50000 

16,5000n 

17 

.80000 

17 

,30000 

17.30000 

17,10000    1 

['.40000 

16, 

90000 

17.J0000 

17 

.60000 

16 

,90000 

16.70000 

16,80000    1 

[0.80000 

17, 

20000 

16,80000 

17, 

60000 

17, 

20000 

16,60000 

17,10000    1 

[6,90000 

16, 

60000 

18,00000 

17, 

20000 

17, 

30000 

17.00000 

16.Q0000    1 

['.30000 

16, 

80000 

17.^0000 

17, 

40000 

17, 

70000 

16,80000 

16,90000    1 

'.OOOOO 

16, 

90000 

17,00000 

16, 

60000 

16, 

70000 

16,80000 

16,70000    1 

[6.AO000 

16, 

50000 

16.40000 

16 

.60000 

16, 

,50000 

16,70000 

16.40000    1 

[6.40000 

16, 

20000 

16.40000 

16 

,30000 

16, 

,40000 

17.00000 

16.90000    1 

['.10000 

17, 

10000 

16, 'OoOo 

16 

-90000 

16, 

,50000 

17.20000 

16,40000    1 

['.00000 

17, 

ooooo 

16,'000o 

16 

'20000 

16 

,60000 

16.90000 

t6. 50000    ! 

16.60000 

16 

,60000 

17,00000 

17 

.10000 

17 

,10000 

16.70000 

16.80000 

[6,30000 

16 

,60000 

16,»0000 

16 

.90000 

17 

.10000 

16.80000 

17.00000 

['.20000 

17 

,30000 

17.  iiOOOo 

17 

'30000 

17, 

.20000 

17.20000 

17.50000    1 

[6,90000 

16, 

,90000 

16. *0000 

17 

>00000 

16 

.50000 

16.70000 

16.80000    1 

16,70000 

16 

,70000 

16.60000 

16 

.50000 

17, 

,00000 

16,70000 

16.70000    ! 

16,90000 

17 

,40000 

17,10000 

17, 

ooooo 

16, 

80000 

17,20000 

17.20000    1 

,'.40000 

17 

,20000 

16. ¥0000 

16, 

•80000 

17, 

ooooo 

17,40000 

17.20000    1 

[',20000 
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